/*
  plot climatology data using gnuplot
  Copyright (C) 2009-2013 Bob Mottram
  bob@sluggish.dyndns.org

  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "gnuplot.h"

gnuplot::gnuplot() {


}

gnuplot::~gnuplot() {

}

int gnuplot::plot_active_stations(
								  string plot_script_filename,
								  string plot_data_filename,
								  string active_stations_image_filename,
								  string equator_dist_image_filename,
								  string title_active_stations,
								  string title_equator_dist,
								  string subtitle, float indent, float vpos,
								  vector<tempdata> &data,
								  vector<long long int> &country_codes,
								  vector<stationdata> &stations,
								  int start_year,
								  int end_year,
								  float min_temp, float max_temp,
								  string population_class,
								  vector<float> &area,
								  int image_width,
								  int image_height)
{
	int max_idx = (end_year - start_year+1);
	vector<long long int> active_stations[max_idx];
	for (int i = 0; i < max_idx; i++) {
		vector<long long int> s;
		active_stations[i] = s;
	}

	int max_active_stations = 1;
	for (int i = 0; i < (int)data.size(); i++) {
		if ((data[i].year >= start_year) &&
			(data[i].year <= end_year) &&
			(data[i].temperature>=min_temp) &&
			(data[i].temperature<=max_temp)) {

			bool is_valid = true;
			if ((int)country_codes.size() > 0) {
                is_valid = ghcn::vector_contains(country_codes, (long long int)data[i].country_code);
			}

			if ((is_valid) && (population_class!="")) {
				is_valid = ghcn::station_has_property(PROPERTY_POPULATION_CLASS,population_class,
													  data[i].station, stations);
			}

			if ((is_valid) && (area.size()==4)) {
				is_valid = ghcn::station_within_area(data[i].station,stations,
													 area[0],area[1],area[2],area[3]);
			}

			if (is_valid) {
			    int idx = data[i].year - start_year;

				if (!ghcn::vector_contains(active_stations[idx], data[i].station)) {
					active_stations[idx].push_back(data[i].station);
				}
			}
		}
	}

	// get distance from the equator
	float dist_from_equator[max_idx];
	float max_dist = 0;
	float min_dist = 90;
    for (int y = start_year; y <= end_year; y++) {
    	int idx = (y - start_year);
    	if ((int)active_stations[idx].size() > max_active_stations) {
    		max_active_stations = (int)active_stations[idx].size();
    	}
    	int hits = 0;
    	dist_from_equator[idx] = 0;
		for (int j = 0; j < (int)active_stations[idx].size(); j++) {
			for (int k = 0; k < (int)stations.size(); k++) {
				if (stations[k].number == active_stations[idx][j]) {
					dist_from_equator[idx] += fabs(stations[k].latitude);
					hits++;
					break;
				}
			}
		}
	    if (hits > 0) {
	    	dist_from_equator[idx] /= hits;
	    	if (dist_from_equator[idx] > max_dist) max_dist = dist_from_equator[idx] + 2;
	    	if (dist_from_equator[idx] < min_dist) min_dist = dist_from_equator[idx] - 2;
	    	if (min_dist < 0) min_dist = 0;
	    }
    }

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    ofstream plotfile;
    plotfile.open(plot_data_filename.c_str());
    for (int y = start_year; y <= end_year; y++) {
    	int idx = y - start_year;
    	plotfile << y << "    " << (int)active_stations[idx].size() << "    " << dist_from_equator[idx] << "\n";
    }
	plotfile.close();

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title_active_stations << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos << "\n";
	plotfile << "set yrange [" << 0 << ":" << (max_active_stations+2) << "]\n";
	plotfile << "set xrange [" << start_year << ":" << end_year << "]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Year\"\n";
	plotfile << "set ylabel \"Number of active stations\"\n";
	plotfile << "set grid\n";
	plotfile << "set key right bottom\n";

	string file_type = "png";
	if (active_stations_image_filename != "") {
		if ((active_stations_image_filename.substr((int)active_stations_image_filename.size()-3,3) == "jpg") ||
			(active_stations_image_filename.substr((int)active_stations_image_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = active_stations_image_filename.substr((int)active_stations_image_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";
		plotfile << "set output \"" << active_stations_image_filename << "\"\n";
		plotfile << "plot \"" << plot_data_filename << "\" using 1:2 notitle with lines\n";
	}

   	if (equator_dist_image_filename != "") {
		if ((equator_dist_image_filename.substr((int)equator_dist_image_filename.size()-3,3) == "jpg") ||
			(equator_dist_image_filename.substr((int)equator_dist_image_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = equator_dist_image_filename.substr((int)equator_dist_image_filename.size()-3,3);
		}

		plotfile << "set title \"" << title_equator_dist << "\"\n";
		plotfile << "set yrange [" << min_dist << ":" << max_dist << "]\n";
		plotfile << "set ylabel \"Average distance to the equator (Degrees)\"\n";
		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";
		plotfile << "set output \"" << equator_dist_image_filename << "\"\n";
		plotfile << "plot \"" << plot_data_filename << "\" using 1:3 notitle with lines\n";
   	}

	plotfile.close();

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}



int gnuplot::plot_annual(
						 string plot_script_filename,
						 string plot_data_filename,
						 string entities_image_filename,
						 string averages_image_filename,
						 string stations_filename,
						 string title,
						 string subtitle, float indent, float vpos,
						 vector<tempdata> &data,
						 int entity_type,
						 vector<long long int> &entity_codes,
						 vector<string> &entity_names,
						 int start_year,
						 int end_year,
						 float min_temp,
						 float max_temp,
						 bool show_minmax,
						 bool best_fit_line,
						 int image_width,
						 int image_height)
{
	const float temperature_margin = 1.0f;

	// keep track of the station codes found
	vector<long long int> station_codes;

	int max_idx = (end_year - start_year+1)*12;
	int no_of_entities = (int)entity_codes.size();
	if ((int)entity_names.size()<no_of_entities) {
		no_of_entities = (int)entity_names.size();
	}

	float data_completeness_percent[no_of_entities+1];
	float annual_max[no_of_entities+1][end_year - start_year+1];
	float annual_min[no_of_entities+1][end_year - start_year+1];
	float month[no_of_entities+1][max_idx];
	int month_hits[no_of_entities+1][max_idx];
	for (int i = 0; i < no_of_entities+1; i++) {
	    for (int j = 0; j < max_idx; j++) {
	    	month[i][j] = -9999;
	    	month_hits[i][j] = 0;
	    }
		for (int j = 0; j < end_year - start_year+1; j++) {
			annual_min[i][j] = 9999;
			annual_max[i][j] = -9999;
		}
		data_completeness_percent[i] = 0;
	}

	float minimum_temp = 9999;
	float maximum_temp = -9999;
	long long int prev_station = -1;
	for (int c = 0; c < no_of_entities+1; c++) {
		for (int i = 0; i < (int)data.size(); i++) {
			int idx = ((data[i].year - start_year)*12) + data[i].month;
			if ((data[i].temperature<min_temp) || (data[i].temperature>max_temp)) continue;
			if (month[no_of_entities][idx] != -9999) {
				month[no_of_entities][idx] += data[i].temperature;
			}
			else {
				month[no_of_entities][idx] = data[i].temperature;
			}
			month_hits[no_of_entities][idx]++;

			if (c != no_of_entities) {
				if (((entity_type == ENTITY_COUNTRY) && (data[i].country_code == entity_codes[c])) ||
					((entity_type == ENTITY_STATION) && (data[i].station == entity_codes[c]))) {
					if (month[c][idx] != -9999) {
						month[c][idx] += data[i].temperature;
					}
					else {
						month[c][idx] = data[i].temperature;
					}
					month_hits[c][idx]++;
					if (data[i].station > 0) {
						if (data[i].station != prev_station) {
							if (!ghcn::vector_contains(station_codes, data[i].station)) {
								station_codes.push_back(data[i].station);
							}
						}
						prev_station = data[i].station;
					}
				}
			}
		}

		int misses = 0;
		for (int i = 0; i < max_idx; i++) {
			if (month_hits[c][i] == 0) {
				misses++;
			}
			else {
				month[c][i] /= month_hits[c][i];
			}
		}
		if ((no_of_entities > 0) && (c < no_of_entities)) {
			if (no_of_entities < 100) {
				if (misses == 0) {
					printf("%s has a complete temperature record\n", entity_names[c].c_str());
				}
				else {
					if ((max_idx > 0) && (max_idx != misses)) {
						printf("%s has a %.1f percent complete temperature record\n", entity_names[c].c_str(),(max_idx-misses)*100.0f/max_idx);
					}
					else {
						printf("%s has no temperature record\n", entity_names[c].c_str());
					}
				}
			}
			if (max_idx > 0) data_completeness_percent[c] = (max_idx-misses)*100.0f/max_idx;
		}
	}

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    for (int y = start_year; y <= end_year; y++) {
    	int idx = (y - start_year) * 12;
	    for (int c = 0; c < no_of_entities+1; c++) {
	    	float annual_average = 0;
	    	int hits = 0;
		    for (int m = 0; m < 12; m++) {

		    	if (month_hits[c][idx+m] > 0) {
		    		float month_temperature = month[c][idx+m];
		    		if (month_temperature > annual_max[c][y-start_year]) {
		    			annual_max[c][y-start_year] = month_temperature;
		    		}
		    		if (month_temperature < annual_min[c][y-start_year]) {
		    			annual_min[c][y-start_year] = month_temperature;
		    		}
		    		annual_average += month_temperature;
		    		hits++;
		    	}
		    }
		    if (hits > 0) {
		    	annual_average /= hits;
		    	if (annual_average < minimum_temp) minimum_temp = annual_average - temperature_margin;
		    	if (annual_average > maximum_temp) maximum_temp = annual_average + temperature_margin;

		    	if (show_minmax) {
			    	if (annual_min[c][y-start_year] < minimum_temp) minimum_temp = annual_min[c][y-start_year] - temperature_margin;
			    	if (annual_max[c][y-start_year] > maximum_temp) maximum_temp = annual_max[c][y-start_year] + temperature_margin;
		    	}
		    }
	    }
    }

    ofstream plotfile;
    plotfile.open(plot_data_filename.c_str());
    for (int y = start_year; y <= end_year; y++) {
    	plotfile << y << "    ";
    	int idx = (y - start_year) * 12;
	    for (int c = 0; c < no_of_entities+1; c++) {
	    	float annual_average = 0;
	    	int hits = 0;
		    for (int m = 0; m < 12; m++) {
		    	if (month_hits[c][idx+m] > 0) {
		    		annual_average += month[c][idx+m];
		    		hits++;
		    	}
		    }
		    if (hits > 0) annual_average /= hits;
		    plotfile << annual_average << "    ";
	    }
	    plotfile << annual_min[no_of_entities][y-start_year] << "    " << annual_max[no_of_entities][y-start_year] << "\n";
    }
	plotfile.close();

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos << "\n";
	plotfile << "set yrange [" << minimum_temp << ":" << maximum_temp << "]\n";
	plotfile << "set xrange [" << start_year << ":" << end_year << "]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Year\"\n";
	plotfile << "set ylabel \"Average Temperature (Celcius)\"\n";
	plotfile << "set grid\n";
	plotfile << "set key right bottom\n";

	string file_type = "png";
	if (entities_image_filename != "") {
		if ((entities_image_filename.substr((int)entities_image_filename.size()-3,3) == "jpg") ||
			(entities_image_filename.substr((int)entities_image_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = entities_image_filename.substr((int)entities_image_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";

		if (no_of_entities > 0) {
			plotfile << "set output \"" << entities_image_filename << "\"\n";
			plotfile << "plot ";

			for (int c = 0; c < no_of_entities; c++) {
				plotfile << "\"" << plot_data_filename << "\" using 1:" << (c+2);
				if ((entity_type == ENTITY_COUNTRY) ||
					((entity_type == ENTITY_STATION) && (data_completeness_percent[c] > 0))) {
					plotfile << " title \"" << entity_names[c] << "\"";
				}
				else {
					plotfile << " notitle";
				}
				plotfile << " with lines";
				if ((c < no_of_entities-1) || (show_minmax)) {
					plotfile << ", ";
				}
				else {
					plotfile << "\n";
				}
			}
			if (show_minmax) {
				plotfile << "\"" << plot_data_filename << "\" using 1:" << (no_of_entities+3) << " title \"Min\" with lines";
				plotfile << ", \"" << plot_data_filename << "\" using 1:" << (no_of_entities+4) << " title \"Max\" with lines\n";
			}
		}
	}

	if (averages_image_filename != "") {
		if ((averages_image_filename.substr((int)averages_image_filename.size()-3,3) == "jpg") ||
			(averages_image_filename.substr((int)averages_image_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = averages_image_filename.substr((int)averages_image_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";

		if (best_fit_line) {
			plotfile << "b = 0.0\n";
			plotfile << "m = 0.0\n";
			plotfile << "f(x) = m*x + b\n";
			plotfile << "fit f(x) \"" << plot_data_filename << "\" using 1:" << (no_of_entities+2) << " via b, m\n";
		}
		plotfile << "set output \"" << averages_image_filename << "\"\n";
		plotfile << "plot \"" << plot_data_filename << "\" using 1:" << (no_of_entities+2);
		if (!show_minmax) {
			plotfile << " notitle with lines";
			if (best_fit_line) plotfile << ", f(x) title \"Best fit\"";
		}
		else {
			plotfile << " title \"Mean\" with lines,";
			plotfile << "\"" << plot_data_filename << "\" using 1:" << (no_of_entities+3) << " title \"Min\" with lines";
			plotfile << ", \"" << plot_data_filename << "\" using 1:" << (no_of_entities+4) << " title \"Max\" with lines";
			if (best_fit_line) {
				plotfile << ", f(x) title \"Best fit\"";
			}
		}
		plotfile << "\n";
	}

	plotfile.close();

	if (stations_filename != "") {
		std::remove(stations_filename.c_str());
	    plotfile.open(stations_filename.c_str());
	    for (int i = 0; i < (int)station_codes.size(); i++) {
	        plotfile << station_codes[i] << "\n";
	    }
	    plotfile.close();
	}

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}


int gnuplot::plot_annual_only(
							  string plot_script_filename,
							  string plot_data_filename,
							  string averages_image_filename,
							  string title,
							  string subtitle, float indent, float vpos,
							  vector<tempdata> &data,
							  vector<long long int> &entity_codes,
							  vector<string> &entity_names,
							  int start_year,
							  int end_year,
							  bool show_minmax,
							  bool best_fit_line,
							  int image_width,
							  int image_height)
{
	// keep track of the station codes found
	vector<long long int> station_codes;

	int max_idx = end_year - start_year+1;
	float average[max_idx];
	int average_hits[max_idx];
	float annual_max[end_year - start_year+1];
	float annual_min[end_year - start_year+1];
	for (int j = 0; j < max_idx; j++) {
		average[j] = -9999;
		average_hits[j] = 0;
	}
	for (int j = 0; j < end_year - start_year+1; j++) {
		annual_min[j] = 9999;
		annual_max[j] = -9999;
	}

	float minimum_temp = 9999;
	float maximum_temp = -9999;
	long long int prev_station = -1;
	for (int i = 0; i < (int)data.size(); i++) {
		int idx = data[i].year - start_year;

		if (average[idx] != -9999) {
			average[idx] += data[i].temperature;
		}
		else {
			average[idx] = data[i].temperature;
		}

		if (data[i].temperature < annual_min[idx]) {
			annual_min[idx] = data[i].temperature;
			if (annual_min[idx] < minimum_temp) minimum_temp = annual_min[idx];
		}
		if (data[i].temperature > annual_max[idx]) {
			annual_max[idx] = data[i].temperature;
			if (annual_max[idx] > maximum_temp) maximum_temp = annual_max[idx];
		}

		average_hits[idx]++;
		if (data[i].station > 0) {
			if (data[i].station != prev_station) {
				if (!ghcn::vector_contains(station_codes, data[i].station)) {
					station_codes.push_back(data[i].station);
				}
			}
			prev_station = data[i].station;
		}
	}

	for (int i = 0; i < max_idx; i++) {
		if (average_hits[i] > 0) {
			average[i] /= average_hits[i];
			if (average[i] < minimum_temp) minimum_temp = average[i];
			if (average[i] > maximum_temp) maximum_temp = average[i];
		}
	}

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    ofstream plotfile;
    plotfile.open(plot_data_filename.c_str());
    for (int y = start_year; y <= end_year; y++) {
    	int idx = y - start_year;
	    plotfile << y << "    " << average[idx] << "    ";
	    plotfile << annual_min[y-start_year] << "    " << annual_max[y-start_year] << "\n";
    }
	plotfile.close();

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos << "\n";
	plotfile << "set yrange [" << minimum_temp << ":" << maximum_temp << "]\n";
	plotfile << "set xrange [" << start_year << ":" << end_year << "]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Year\"\n";
	plotfile << "set ylabel \"Average Temperature (Celcius)\"\n";
	plotfile << "set grid\n";
	plotfile << "set key right bottom\n";

	string file_type = "png";
	if (averages_image_filename != "") {
		if ((averages_image_filename.substr((int)averages_image_filename.size()-3,3) == "jpg") ||
			(averages_image_filename.substr((int)averages_image_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = averages_image_filename.substr((int)averages_image_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";

		if (best_fit_line) {
			plotfile << "b = 0.0\n";
			plotfile << "m = 0.0\n";
			plotfile << "f(x) = m*x + b\n";
			plotfile << "fit f(x) \"" << plot_data_filename << "\" using 1:2 via b, m\n";
		}
		plotfile << "set output \"" << averages_image_filename << "\"\n";
		plotfile << "plot \"" << plot_data_filename << "\" using 1:2";
		if (!show_minmax) {
			plotfile << " notitle with lines";
			if (best_fit_line) plotfile << ", f(x) title \"Best fit\"";
		}
		else {
			plotfile << " title \"Mean\" with lines,";
			plotfile << "\"" << plot_data_filename << "\" using 1:3 title \"Min\" with lines";
			plotfile << ", \"" << plot_data_filename << "\" using 1:4 title \"Max\" with lines";
			if (best_fit_line) {
				plotfile << ", f(x) title \"Best fit\"";
			}
		}
		plotfile << "\n";
	}

	plotfile.close();

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}


void gnuplot::plot_by_country_annual(
									 string plot_script_filename,
									 string plot_data_filename,
									 string countries_image_filename,
									 string averages_image_filename,
									 string stations_filename,
									 string title,
									 string subtitle, float indent, float vpos,
									 vector<tempdata> &data,
									 vector<long long int> &country_codes,
									 vector<string> &country_names,
									 int start_year,
									 int end_year,
									 float min_temp,
									 float max_temp,
									 bool show_minmax,
									 bool best_fit_line,
									 int image_width,
									 int image_height)
{
	plot_annual(
				plot_script_filename,
				plot_data_filename,
				countries_image_filename,
				averages_image_filename,
				stations_filename,
				title,
				subtitle, indent, vpos,
				data,
				ENTITY_COUNTRY,
				country_codes,
				country_names,
				start_year,
				end_year,
				min_temp,
				max_temp,
				show_minmax,
				best_fit_line,
				image_width,
				image_height);
}


void gnuplot::plot_by_station_annual(
									 string plot_script_filename,
									 string plot_data_filename,
									 string stations_image_filename,
									 string averages_image_filename,
									 string title,
									 string subtitle, float indent, float vpos,
									 vector<tempdata> &data,
									 vector<long long int> &station_numbers,
									 vector<string> &station_names,
									 int start_year,
									 int end_year,
									 bool show_minmax,
									 bool best_fit_line,
									 int image_width,
									 int image_height)
{
	plot_annual(
				plot_script_filename,
				plot_data_filename,
				stations_image_filename,
				averages_image_filename,
				"",
				title,
				subtitle, indent, vpos,
				data,
				ENTITY_STATION,
				station_numbers,
				station_names,
				start_year,
				end_year,
				-100.0f, 100.0f,
				show_minmax,
				best_fit_line,
				image_width,
				image_height);
}


int gnuplot::plot_anomalies(
							string plot_script_filename,
							string plot_data_filename,
							string anomalies_image_filename,
							string title,
							string subtitle, float indent, float vpos,
							vector<gridcell> &grid,
							int start_year,
							int end_year,
							int reference_start_year,
							int reference_end_year,
							bool show_minmax,
							bool show_running_average,
							bool best_fit_line,
							int image_width,
							int image_height)
{
	float minimum_anomaly=9999,maximum_anomaly=-9999;
	vector<float> time_series;
	vector<float> time_series_min;
	vector<float> time_series_max;

	printf("Calculating reference temperatures...");
	
	for (int i = 0; i < (int)grid.size(); i++) {
		grid[i].update(reference_start_year,reference_end_year);
	}

	printf("done\nPlotting anomalies...");

	for (int year=start_year; year<=end_year; year++) {
		float anomaly=0,av_min=0,av_max=0;
		int anomaly_hits=0,av_hits=0;
		for (int i = 0; i < (int)grid.size(); i++) {
			float anom = grid[i].get_anomaly(year,reference_start_year,reference_end_year);
			if (anom>-9999.0f) {
				float min=9999,max=-9999;
				if (show_minmax) {
					grid[i].get_variance(year,min,max,reference_start_year,reference_end_year);
					av_min += min;
					av_max += max;
					av_hits++;
				}
				anomaly += anom;
				anomaly_hits++;
			}
		}
		if (anomaly_hits>0) {
			anomaly /= anomaly_hits;
			if (show_minmax) {
				av_min = (float)(av_min / av_hits);
				av_max = (float)(av_max / av_hits);
				if (av_min < minimum_anomaly) minimum_anomaly = av_min;
				if (av_max > maximum_anomaly) maximum_anomaly = av_max;
			}
			else {
				if (anomaly < minimum_anomaly) minimum_anomaly = anomaly;
				if (anomaly > maximum_anomaly) maximum_anomaly = anomaly;
			}
		}
		time_series.push_back(anomaly);
		time_series_min.push_back(av_min);
		time_series_max.push_back(av_max);
	}

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    ofstream plotfile;
    plotfile.open(plot_data_filename.c_str());
	float running_average = time_series[0];
	if ((show_running_average) && (!show_minmax)) {
		minimum_anomaly=9999;
		maximum_anomaly=-9999;
	}
    for (int y = start_year; y <= end_year; y++) {
    	int idx = y - start_year;
		plotfile << y << "    " << time_series[idx];
		plotfile << "    " << time_series_min[idx];
		plotfile << "    " << time_series_max[idx];
		plotfile << "    " << running_average;
		plotfile << "\n";
		running_average = (running_average*0.8f)+(time_series[idx]*0.2f);

		if ((show_running_average) && (!show_minmax)) {
			if (running_average>maximum_anomaly) maximum_anomaly = running_average;
			if (running_average<minimum_anomaly) minimum_anomaly = running_average;
		}
    }
	plotfile.close();

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos <<"\n";
	plotfile << "set yrange [" << minimum_anomaly << ":" << maximum_anomaly << "]\n";
	plotfile << "set xrange [" << start_year << ":" << end_year << "]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Year\"\n";
	plotfile << "set ylabel \"Average Anomaly (Celcius) relative to ";
	plotfile << reference_start_year << " - " << reference_end_year << "\"\n";
	plotfile << "set grid\n";
	plotfile << "set key right bottom\n";

	string file_type = "png";
	if (anomalies_image_filename != "") {
		if ((anomalies_image_filename.substr((int)anomalies_image_filename.size()-3,3) == "jpg") ||
			(anomalies_image_filename.substr((int)anomalies_image_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = anomalies_image_filename.substr((int)anomalies_image_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";

		if (best_fit_line) {
			plotfile << "b = 0.0\n";
			plotfile << "m = 0.0\n";
			plotfile << "f(x) = m*x + b\n";
			plotfile << "fit f(x) \"" << plot_data_filename << "\" using 1:2 via b, m\n";
		}
		plotfile << "set output \"" << anomalies_image_filename << "\"\n";
		if (show_running_average) {
			plotfile << "plot \"" << plot_data_filename << "\" using 1:5";
		}
		else {
			plotfile << "plot \"" << plot_data_filename << "\" using 1:2";
		}

		if (!show_minmax) {
			plotfile << " notitle with lines";
			if (best_fit_line) plotfile << ", f(x) title \"Best fit\"";
		}
		else {
			plotfile << " title \"Mean\" with lines,";
			plotfile << "\"" << plot_data_filename << "\" using 1:3 title \"Min\" with lines";
			plotfile << ", \"" << plot_data_filename << "\" using 1:4 title \"Max\" with lines";
			if (best_fit_line) {
				plotfile << ", f(x) title \"Best fit\"";
			}
		}


		plotfile << "\n";
	}

	plotfile.close();

	printf("done\n");

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}


int gnuplot::plot_distribution(
							   string plot_script_filename,
							   string plot_data_filename,
							   string anomalies_image_filename,
							   string title,
							   string subtitle, float indent, float vpos,
							   vector<gridcell> &grid,
							   int start_year,
							   int end_year,
							   int reference_start_year,
							   int reference_end_year,
							   int interval,
							   int image_width,
							   int image_height)
{
	int intervals = 1;
	if (interval>0) intervals = (end_year-start_year)/interval;
	if (intervals<1) intervals=1;
	float histogram[intervals][200];
	int hits[intervals];
	int min_temp=9999,max_temp=-9999;

	printf("Plotting distribution...");

	for (int p = 0; p < intervals; p++) {

		hits[p]=0;
		for (int s=0; s < 200; s++) histogram[p][s]=0;

		int interval_start_year = start_year + (p*(end_year-start_year)/intervals);
		int interval_end_year = start_year + ((p+1)*(end_year-start_year)/intervals);

		for (int year=interval_start_year; year<interval_end_year; year++) {
			for (int i = 0; i < (int)grid.size(); i++) {
				if (grid[i].update_histogram(year)) {
					for (int s = 0; s < 200; s++) {
						histogram[p][s] += grid[i].histogram[s];
					}
					hits[p]++;
				}
			}
		}

		if (hits[p]>0) {
			for (int i = 0; i < 200; i++) {
				histogram[p][i] /= (float)hits[p];
				if (histogram[p][i]>0) {
					if (i-100>max_temp) max_temp = i-100;
					if (i-100<min_temp) min_temp = i-100;
				}		
			}
		}
	}

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    ofstream plotfile;
    plotfile.open(plot_data_filename.c_str());
	float max=0;
    for (int s = 0; s < 200; s++) {
		plotfile << (s-100);
		for (int p = 0; p < intervals; p++) {
			plotfile << "    " << histogram[p][s];
			if (histogram[p][s]>max) max = histogram[p][s];
		}
		plotfile << "\n";
    }
	plotfile.close();

	if (max==0) {
		printf("No data\n");
		return -1;
	}

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos << "\n";
	plotfile << "set yrange [0:" << max << "]\n";
	plotfile << "set xrange [" << min_temp << ":" << max_temp << "]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Temperature (Celcius)\"\n";
	plotfile << "set ylabel \"Normalised Samples\"\n";
	plotfile << "set grid\n";
	plotfile << "set key left top\n";

	string file_type = "png";
	if (anomalies_image_filename != "") {
		if ((anomalies_image_filename.substr((int)anomalies_image_filename.size()-3,3) == "jpg") ||
			(anomalies_image_filename.substr((int)anomalies_image_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = anomalies_image_filename.substr((int)anomalies_image_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";

		plotfile << "set output \"" << anomalies_image_filename << "\"\n";

		plotfile << "plot ";

		for (int p = 0; p < intervals; p++) {
			if (p>0) plotfile << ",";
			plotfile << "\"" << plot_data_filename << "\" using 1:" << (2+p);
			if (intervals>1) {
				plotfile << " title \"" << (start_year+(p*interval)) << " - " << (start_year+((p+1)*interval)) << "\"";
			}
			else {
				plotfile << " notitle";
			}
			plotfile << " with lines";
		}

		plotfile << "\n";
	}
	else {
		printf("\nERROR: No image filename specified\n");
	}

	plotfile.close();

	printf("done\n");

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}


int gnuplot::plot_seasonal_anomalies(
									 string plot_script_filename,
									 string plot_data_filename,
									 string anomalies_image_filename,
									 string title,
									 string subtitle, float indent, float vpos,
									 vector<gridcell> &grid,
									 int start_year,
									 int end_year,
									 int reference_start_year,
									 int reference_end_year,
									 bool show_running_average,
									 int image_width,
									 int image_height)
{
	float minimum_anomaly=9999,maximum_anomaly=-9999;
	vector<float> time_series[12];

	printf("Calculating reference temperatures...");
	
	for (int i = 0; i < (int)grid.size(); i++) {
		grid[i].update(reference_start_year,reference_end_year);
	}

	printf("done\nPlotting seasonal anomalies...");

	for (int year=start_year; year<=end_year; year++) {		
		for (int m = 0; m < 12; m++) {
			float anomaly=0;
			int anomaly_hits=0;
			for (int i = 0; i < (int)grid.size(); i++) {
				float anom = grid[i].get_anomaly(year,m,reference_start_year,reference_end_year);
				if (anom>-9999.0f) {
					anomaly += anom;
					anomaly_hits++;
				}
			}
		
			if (anomaly_hits>0) {
				anomaly /= anomaly_hits;
			}
			time_series[m].push_back(anomaly);
		}
	}

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    ofstream plotfile;
    plotfile.open(plot_data_filename.c_str());
	float running_average[2];
	running_average[0] = (time_series[11][0]+time_series[0][0]+time_series[1][0]+time_series[2][0])/4;
	running_average[1] = (time_series[5][0]+time_series[6][0]+time_series[7][0]+time_series[8][0])/4;

    for (int y = start_year; y <= end_year; y++) {
    	int idx = y - start_year;
		plotfile << y;

		// winter Dec,Jan,Feb
		float winter = (time_series[11][idx]+time_series[0][idx]+time_series[1][idx]+time_series[2][idx])/4;
		plotfile << "    " << winter;
		plotfile << "    " << running_average[0];
		running_average[0] =
				(running_average[0]*0.8f)+(winter*0.2f);

		// Summer Jun, Jul, Aug
		float summer = (time_series[5][idx]+time_series[6][idx]+time_series[7][idx]+time_series[8][idx])/4;
		plotfile << "    " << summer;
		plotfile << "    " << running_average[1];
		running_average[1] =
				(running_average[1]*0.8f)+(summer*0.2f);

		plotfile << "\n";

		if (show_running_average) {
			if (running_average[0] < minimum_anomaly) minimum_anomaly = running_average[0];
			if (running_average[0] > maximum_anomaly) maximum_anomaly = running_average[0];
			if (running_average[1] < minimum_anomaly) minimum_anomaly = running_average[1];
			if (running_average[1] > maximum_anomaly) maximum_anomaly = running_average[1];
		}
		else {
			if (winter < minimum_anomaly) minimum_anomaly = winter;
			if (winter > maximum_anomaly) maximum_anomaly = winter;
			if (summer < minimum_anomaly) minimum_anomaly = summer;
			if (summer > maximum_anomaly) maximum_anomaly = summer;
		}
    }
	plotfile.close();

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos << "\n";
	plotfile << "set yrange [" << minimum_anomaly << ":" << maximum_anomaly << "]\n";
	plotfile << "set xrange [" << start_year << ":" << end_year << "]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Year\"\n";
	plotfile << "set ylabel \"Average Anomaly (Celcius) relative to ";
	plotfile << reference_start_year << " - " << reference_end_year << "\"\n";
	plotfile << "set grid\n";
	plotfile << "set key right bottom\n";

	string file_type = "png";
	if (anomalies_image_filename != "") {
		if ((anomalies_image_filename.substr((int)anomalies_image_filename.size()-3,3) == "jpg") ||
			(anomalies_image_filename.substr((int)anomalies_image_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = anomalies_image_filename.substr((int)anomalies_image_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";

		plotfile << "set output \"" << anomalies_image_filename << "\"\n";
		plotfile << "plot ";

		if (show_running_average) {
			plotfile << "\"" << plot_data_filename << "\" using 1:3 title \"Dec-Mar\" with lines,";
			plotfile << "\"" << plot_data_filename << "\" using 1:5 title \"Jun-Sep\" with lines";
		}
		else {
			plotfile << "\"" << plot_data_filename << "\" using 1:2 title \"Dec-Mar\" with lines,";
			plotfile << "\"" << plot_data_filename << "\" using 1:4 title \"Jun-Sep\" with lines";
		}
		plotfile << "\n";
	}

	plotfile.close();

	printf("done\n");

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}

void gnuplot::save_world_map(string filename)
{
    ofstream plotfile;

	std::remove(filename.c_str());

    plotfile.open(filename.c_str());

	plotfile << "#\n";
	plotfile << "# $Id: world.dat,v 1.2 2006/06/02 06:01:44 sfeam Exp $\n";
	plotfile << "#\n";
	plotfile << "#\n";
	plotfile << "-92.32  48.24\n";
	plotfile << "-88.13  48.92\n";
	plotfile << "-83.11  46.27\n";
	plotfile << "-81.66  44.76\n";
	plotfile << "-82.09  42.29\n";
	plotfile << "-77.10  44.00\n";
	plotfile << "-69.95  46.92\n";
	plotfile << "-65.92  45.32\n";
	plotfile << "-66.37  44.25\n";
	plotfile << "-61.22  45.43\n";
	plotfile << "-64.94  47.34\n";
	plotfile << "-64.12  48.52\n";
	plotfile << "-70.68  47.02\n";
	plotfile << "-67.24  49.33\n";
	plotfile << "-59.82  50.48\n";
	plotfile << "-56.14  52.46\n";
	plotfile << "-59.07  53.58\n";
	plotfile << "-58.26  54.21\n";
	plotfile << "-60.69  55.33\n";
	plotfile << "-61.97  57.41\n";
	plotfile << "-64.35  59.49\n";
	plotfile << "-67.29  58.15\n";
	plotfile << "-69.89  59.91\n";
	plotfile << "-71.31  61.45\n";
	plotfile << "-78.22  61.97\n";
	plotfile << "-77.28  59.53\n";
	plotfile << "-77.09  55.88\n";
	plotfile << "-79.06  51.68\n";
	plotfile << "-82.23  52.70\n";
	plotfile << "-86.75  55.72\n";
	plotfile << "-92.17  56.86\n";
	plotfile << "-95.61  58.82\n";
	plotfile << "-92.66  62.02\n";
	plotfile << "-90.65  63.24\n";
	plotfile << "-95.96  64.12\n";
	plotfile << "-89.88  63.98\n";
	plotfile << "-89.30  65.22\n";
	plotfile << "-86.86  66.12\n";
	plotfile << "-84.54  66.88\n";
	plotfile << "-82.30  67.76\n";
	plotfile << "-83.10  69.68\n";
	plotfile << "-86.05  67.98\n";
	plotfile << "-88.18  68.20\n";
	plotfile << "-91.00  68.82\n";
	plotfile << "-91.72  69.69\n";
	plotfile << "-93.15  71.09\n";
	plotfile << "-96.58  71.05\n";
	plotfile << "-93.35  69.52\n";
	plotfile << "-94.23  68.25\n";
	plotfile << "-95.96  66.73\n";
	plotfile << "-98.83  68.27\n";
	plotfile << "-102.45  67.69\n";
	plotfile << "-108.34  68.43\n";
	plotfile << "-105.83  68.05\n";
	plotfile << "-108.15  66.60\n";
	plotfile << "-111.15  67.63\n";
	plotfile << "-114.10  68.23\n";
	plotfile << "-120.92  69.44\n";
	plotfile << "-124.32  69.26\n";
	plotfile << "-128.76  70.50\n";
	plotfile << "-131.86  69.19\n";
	plotfile << "-131.15  69.79\n";
	plotfile << "-135.81  69.13\n";
	plotfile << "-140.19  69.37\n";
	plotfile << "-141.20  69.58 \n";
	plotfile << "-141.21  69.56 \n";
	plotfile << "-142.49  69.83\n";
	plotfile << "-148.09  70.26\n";
	plotfile << "-154.37  70.96\n";
	plotfile << "-159.53  70.38\n";
	plotfile << "-166.64  68.25\n";
	plotfile << "-161.56  66.55\n";
	plotfile << "-162.99  65.97\n";
	plotfile << "-168.23  65.49\n";
	plotfile << "-161.12  64.49\n";
	plotfile << "-165.29  62.57\n";
	plotfile << "-164.58  60.06\n";
	plotfile << "-162.06  58.36\n";
	plotfile << "-157.85  58.12\n";
	plotfile << "-162.34  55.06\n";
	plotfile << "-156.52  57.11\n";
	plotfile << "-153.53  59.32\n";
	plotfile << "-149.18  60.81\n";
	plotfile << "-149.90  59.50\n";
	plotfile << "-146.54  60.36\n";
	plotfile << "-139.98  59.73\n";
	plotfile << "-137.12  58.28\n";
	plotfile << "-136.01  59.12\n";
	plotfile << "-133.84  57.12\n";
	plotfile << "-131.46  55.98\n";
	plotfile << "-132.08  57.20\n";
	plotfile << "-140.37  60.25\n";
	plotfile << "-141.21  60.16\n";
	plotfile << "-133.38  58.93\n";
	plotfile << "-130.88  54.83\n";
	plotfile << "-128.86  53.90\n";
	plotfile << "-126.58  52.12\n";
	plotfile << "-127.08  50.80\n";
	plotfile << "-124.42  49.66\n";
	plotfile << "-122.56  48.91 \n";
	plotfile << "-122.44  48.92\n";
	plotfile << "-124.42  47.18\n";
	plotfile << "-124.52  42.48\n";
	plotfile << "-123.09  38.45\n";
	plotfile << "-121.73  36.62\n";
	plotfile << "-117.60  33.34\n";
	plotfile << "-117.28  32.64 \n";
	plotfile << "-117.29  32.48\n";
	plotfile << "-114.75  27.80\n";
	plotfile << "-112.53  24.80\n";
	plotfile << "-110.55  24.07\n";
	plotfile << "-114.23  29.59\n";
	plotfile << "-112.58  29.99\n";
	plotfile << "-109.57  25.94\n";
	plotfile << "-105.61  21.94\n";
	plotfile << "-102.09  17.87\n";
	plotfile << "-95.75  15.94\n";
	plotfile << "-92.21  14.97 \n";
	plotfile << "-92.22  14.71\n";
	plotfile << "-86.74  12.06\n";
	plotfile << "-83.03   8.65\n";
	plotfile << "-79.93   8.74\n";
	plotfile << "-77.00   7.82\n";
	plotfile << "-81.99   8.97\n";
	plotfile << "-83.92  12.70\n";
	plotfile << "-86.33  15.80\n";
	plotfile << "-88.40  15.92 \n";
	plotfile << "-88.45  17.42\n";
	plotfile << "-87.01  21.33\n";
	plotfile << "-91.65  18.72\n";
	plotfile << "-96.96  20.37\n";
	plotfile << "-97.65  25.67 \n";
	plotfile << "-97.62  25.82\n";
	plotfile << "-95.62  28.84\n";
	plotfile << "-90.77  29.03\n";
	plotfile << "-87.33  30.22\n";
	plotfile << "-82.69  28.15\n";
	plotfile << "-80.16  26.66\n";
	plotfile << "-80.74  32.31\n";
	plotfile << "-76.89  35.43\n";
	plotfile << "-76.47  38.21\n";
	plotfile << "-75.66  37.67\n";
	plotfile << "-71.31  41.76\n";
	plotfile << "-69.44  44.17\n";
	plotfile << "-67.69  47.03\n";
	plotfile << "-73.18  45.14\n";
	plotfile << "-79.26  43.28\n";
	plotfile << "-82.84  42.59\n";
	plotfile << "-83.49  45.32\n";
	plotfile << "-86.36  43.65\n";
	plotfile << "-87.75  43.42\n";
	plotfile << "-86.01  45.96\n";
	plotfile << "-87.00  46.59\n";
	plotfile << "-91.39  46.79\n";
	plotfile << "-90.05  47.96 \n";
	plotfile << "\n";
	plotfile << "-152.62  58.41\n";
	plotfile << "-152.60  58.40 \n";
	plotfile << "\n";
	plotfile << "-153.30  57.80\n";
	plotfile << "-152.40  57.48\n";
	plotfile << "-153.32  57.79 \n";
	plotfile << "\n";
	plotfile << "-166.96  53.96\n";
	plotfile << "-167.01  53.95 \n";
	plotfile << "\n";
	plotfile << "-168.36  53.50\n";
	plotfile << "-168.19  53.36 \n";
	plotfile << "\n";
	plotfile << "-170.73  52.68\n";
	plotfile << "-170.60  52.55 \n";
	plotfile << "\n";
	plotfile << "-174.47  51.94\n";
	plotfile << "-174.47  51.92 \n";
	plotfile << "\n";
	plotfile << "-176.58  51.71\n";
	plotfile << "-176.64  51.73 \n";
	plotfile << "\n";
	plotfile << "-177.55  51.76\n";
	plotfile << "-177.41  51.63 \n";
	plotfile << "\n";
	plotfile << "-178.27  51.75 \n";
	plotfile << "\n";
	plotfile << "177.35  51.80\n";
	plotfile << "177.33  51.76 \n";
	plotfile << "\n";
	plotfile << "172.44  53.00\n";
	plotfile << "172.55  53.03 \n";
	plotfile << "\n";
	plotfile << "-123.40  48.33\n";
	plotfile << "-128.00  50.84\n";
	plotfile << "-123.50  48.34 \n";
	plotfile << "\n";
	plotfile << "-132.49  52.88\n";
	plotfile << "-132.44  52.91 \n";
	plotfile << "\n";
	plotfile << "-132.64  53.02\n";
	plotfile << "-131.97  53.71\n";
	plotfile << "-132.63  53.02 \n";
	plotfile << "\n";
	plotfile << "-55.36  51.56\n";
	plotfile << "-54.66  49.52\n";
	plotfile << "-53.65  47.48\n";
	plotfile << "-52.98  46.31\n";
	plotfile << "-56.12  46.84\n";
	plotfile << "-58.47  47.57\n";
	plotfile << "-57.61  50.38\n";
	plotfile << "-55.39  51.53 \n";
	plotfile << "\n";
	plotfile << "-61.37  49.01\n";
	plotfile << "-61.80  49.29\n";
	plotfile << "-61.38  49.03 \n";
	plotfile << "\n";
	plotfile << "-63.01  46.71\n";
	plotfile << "-64.42  46.61\n";
	plotfile << "-63.04  46.68 \n";
	plotfile << "\n";
	plotfile << "-60.14  46.48\n";
	plotfile << "-60.14  46.50 \n";
	plotfile << "\n";
	plotfile << "-71.97  41.11\n";
	plotfile << "-71.97  41.15 \n";
	plotfile << "\n";
	plotfile << "-80.79  27.03\n";
	plotfile << "-81.01  26.99 \n";
	plotfile << "\n";
	plotfile << "-113.01  42.09\n";
	plotfile << "-113.10  42.01 \n";
	plotfile << "\n";
	plotfile << "-155.74  20.02\n";
	plotfile << "-155.73  19.98 \n";
	plotfile << "\n";
	plotfile << "-156.51  20.78\n";
	plotfile << "-156.51  20.78 \n";
	plotfile << "\n";
	plotfile << "-157.12  21.21\n";
	plotfile << "-157.08  20.95 \n";
	plotfile << "\n";
	plotfile << "-157.87  21.42 \n";
	plotfile << "\n";
	plotfile << "-159.53  22.07 \n";
	plotfile << "\n";
	plotfile << "-117.44  66.46\n";
	plotfile << "-119.59  65.24\n";
	plotfile << "-123.95  65.03\n";
	plotfile << "-123.69  66.44\n";
	plotfile << "-119.21  66.22\n";
	plotfile << "-117.44  66.44 \n";
	plotfile << "\n";
	plotfile << "-120.71  64.03\n";
	plotfile << "-114.91  62.30\n";
	plotfile << "-109.07  62.72\n";
	plotfile << "-112.62  61.19\n";
	plotfile << "-118.68  61.19\n";
	plotfile << "-117.01  61.17\n";
	plotfile << "-115.97  62.56\n";
	plotfile << "-119.46  64.00\n";
	plotfile << "-120.59  63.94 \n";
	plotfile << "\n";
	plotfile << "-112.31  58.46\n";
	plotfile << "-108.90  59.44\n";
	plotfile << "-104.14  58.90\n";
	plotfile << "-102.56  56.72\n";
	plotfile << "-101.82  58.73\n";
	plotfile << "-104.65  58.91\n";
	plotfile << "-111.00  58.51\n";
	plotfile << "-112.35  58.62 \n";
	plotfile << "\n";
	plotfile << "-98.74  50.09\n";
	plotfile << "-99.75  52.24\n";
	plotfile << "-99.62  51.47\n";
	plotfile << "-98.82  50.39 \n";
	plotfile << "\n";
	plotfile << "-97.02  50.21\n";
	plotfile << "-97.50  54.02\n";
	plotfile << "-98.69  52.93\n";
	plotfile << "-97.19  51.09\n";
	plotfile << "-96.98  50.20 \n";
	plotfile << "\n";
	plotfile << "-95.34  49.04\n";
	plotfile << "-92.32  50.34\n";
	plotfile << "-94.14  49.47\n";
	plotfile << "-95.36  48.82 \n";
	plotfile << "\n";
	plotfile << "-80.39  56.16\n";
	plotfile << "-79.22  55.94\n";
	plotfile << "-80.34  56.08 \n";
	plotfile << "\n";
	plotfile << "-103.56  58.60\n";
	plotfile << "-103.60  58.58 \n";
	plotfile << "\n";
	plotfile << "-101.82  58.03\n";
	plotfile << "-102.33  58.10\n";
	plotfile << "-101.77  58.06 \n";
	plotfile << "\n";
	plotfile << "-101.88  55.79\n";
	plotfile << "-97.92  57.15\n";
	plotfile << "-101.22  55.85\n";
	plotfile << "-101.88  55.74 \n";
	plotfile << "\n";
	plotfile << "-77.61   6.80\n";
	plotfile << "-78.70   0.97\n";
	plotfile << "-80.75  -4.47\n";
	plotfile << "-76.19 -14.57\n";
	plotfile << "-70.44 -18.75\n";
	plotfile << "-70.68 -26.15\n";
	plotfile << "-71.44 -32.03\n";
	plotfile << "-73.38 -37.27\n";
	plotfile << "-73.06 -42.11\n";
	plotfile << "-73.17 -46.09\n";
	plotfile << "-73.52 -48.05\n";
	plotfile << "-73.67 -51.56\n";
	plotfile << "-71.06 -53.88\n";
	plotfile << "-69.14 -50.77\n";
	plotfile << "-67.51 -46.59\n";
	plotfile << "-63.49 -42.80\n";
	plotfile << "-62.14 -40.16\n";
	plotfile << "-57.12 -36.71\n";
	plotfile << "-53.17 -34.15\n";
	plotfile << "-51.26 -32.02\n";
	plotfile << "-48.16 -25.48\n";
	plotfile << "-40.73 -22.32\n";
	plotfile << "-38.88 -15.24\n";
	plotfile << "-34.60  -7.81\n";
	plotfile << "-41.95  -3.42\n";
	plotfile << "-48.02  -1.84\n";
	plotfile << "-48.44  -1.57\n";
	plotfile << "-50.81   0.00\n";
	plotfile << "-54.47   5.39\n";
	plotfile << "-60.59   8.32\n";
	plotfile << "-64.19   9.88\n";
	plotfile << "-70.78  10.64\n";
	plotfile << "-70.97  11.89\n";
	plotfile << "-76.26   8.76\n";
	plotfile << "-77.61   6.80 \n";
	plotfile << "\n";
	plotfile << "-69.14 -52.79\n";
	plotfile << "-66.16 -55.08\n";
	plotfile << "-70.01 -54.88\n";
	plotfile << "-70.55 -53.85\n";
	plotfile << "-69.31 -52.81 \n";
	plotfile << "\n";
	plotfile << "-59.29 -51.58\n";
	plotfile << "-59.35 -51.54 \n";
	plotfile << "\n";
	plotfile << "-58.65 -51.55\n";
	plotfile << "-58.55 -51.56 \n";
	plotfile << "\n";
	plotfile << "-84.39  21.44\n";
	plotfile << "-73.90  19.73\n";
	plotfile << "-79.27  21.18\n";
	plotfile << "-83.74  21.80\n";
	plotfile << "-84.32  21.42 \n";
	plotfile << "\n";
	plotfile << "-66.96  17.95\n";
	plotfile << "-67.05  17.89 \n";
	plotfile << "\n";
	plotfile << "-77.88  17.22\n";
	plotfile << "-78.06  16.98 \n";
	plotfile << "\n";
	plotfile << "-74.47  18.08\n";
	plotfile << "-69.88  18.99\n";
	plotfile << "-71.10  17.76\n";
	plotfile << "-74.45  17.86 \n";
	plotfile << "\n";
	plotfile << "-85.28  73.74\n";
	plotfile << "-85.79  70.96\n";
	plotfile << "-85.13  71.94\n";
	plotfile << "-84.74  72.96\n";
	plotfile << "-80.61  73.10\n";
	plotfile << "-78.45  72.20\n";
	plotfile << "-75.44  72.55\n";
	plotfile << "-73.89  71.98\n";
	plotfile << "-72.56  71.04\n";
	plotfile << "-71.49  70.57\n";
	plotfile << "-69.78  70.29\n";
	plotfile << "-68.12  69.71\n";
	plotfile << "-65.91  69.19\n";
	plotfile << "-66.92  68.39\n";
	plotfile << "-64.08  67.68\n";
	plotfile << "-62.50  66.68\n";
	plotfile << "-63.07  65.33\n";
	plotfile << "-66.11  66.08\n";
	plotfile << "-67.48  65.41\n";
	plotfile << "-64.05  63.15\n";
	plotfile << "-66.58  63.26\n";
	plotfile << "-69.04  62.33\n";
	plotfile << "-72.22  63.77\n";
	plotfile << "-76.88  64.17\n";
	plotfile << "-73.25  65.54\n";
	plotfile << "-70.09  66.64\n";
	plotfile << "-72.05  67.44\n";
	plotfile << "-76.32  68.36\n";
	plotfile << "-78.34  70.17\n";
	plotfile << "-82.12  69.71\n";
	plotfile << "-87.64  70.12\n";
	plotfile << "-89.68  71.43\n";
	plotfile << "-85.28  73.74 \n";
	plotfile << "\n";
	plotfile << "-80.90  76.10\n";
	plotfile << "-84.21  76.28\n";
	plotfile << "-88.94  76.38\n";
	plotfile << "-85.47  77.40\n";
	plotfile << "-85.43  77.93\n";
	plotfile << "-87.01  78.54\n";
	plotfile << "-83.17  78.94\n";
	plotfile << "-84.87  79.93\n";
	plotfile << "-81.33  79.82\n";
	plotfile << "-76.27  80.92\n";
	plotfile << "-82.88  80.62\n";
	plotfile << "-82.58  81.16\n";
	plotfile << "-86.51  81.05\n";
	plotfile << "-89.36  81.21\n";
	plotfile << "-90.45  81.38\n";
	plotfile << "-89.28  81.86\n";
	plotfile << "-87.21  82.30\n";
	plotfile << "-80.51  82.05\n";
	plotfile << "-80.16  82.55\n";
	plotfile << "-77.83  82.86\n";
	plotfile << "-75.51  83.05\n";
	plotfile << "-71.18  82.90\n";
	plotfile << "-65.10  82.78\n";
	plotfile << "-63.34  81.80\n";
	plotfile << "-68.26  81.26\n";
	plotfile << "-69.46  80.34\n";
	plotfile << "-71.05  79.82\n";
	plotfile << "-74.40  79.46\n";
	plotfile << "-75.42  79.03\n";
	plotfile << "-75.48  78.92\n";
	plotfile << "-76.01  78.20\n";
	plotfile << "-80.66  77.28\n";
	plotfile << "-78.07  76.98\n";
	plotfile << "-80.90  76.13 \n";
	plotfile << "\n";
	plotfile << "-92.86  74.13\n";
	plotfile << "-92.50  72.70\n";
	plotfile << "-94.89  73.16\n";
	plotfile << "-92.96  74.14 \n";
	plotfile << "\n";
	plotfile << "-94.80  76.95\n";
	plotfile << "-89.68  76.04\n";
	plotfile << "-88.52  75.40\n";
	plotfile << "-82.36  75.67\n";
	plotfile << "-79.39  74.65\n";
	plotfile << "-86.15  74.22\n";
	plotfile << "-91.70  74.94\n";
	plotfile << "-95.60  76.91\n";
	plotfile << "-94.87  76.96 \n";
	plotfile << "\n";
	plotfile << "-99.96  73.74\n";
	plotfile << "-97.89  72.90\n";
	plotfile << "-98.28  71.13\n";
	plotfile << "-102.04  72.92\n";
	plotfile << "-101.34  73.14\n";
	plotfile << "-99.69  73.59 \n";
	plotfile << "\n";
	plotfile << "-107.58  73.25\n";
	plotfile << "-104.59  71.02\n";
	plotfile << "-101.71  69.56\n";
	plotfile << "-104.07  68.62\n";
	plotfile << "-106.61  69.12\n";
	plotfile << "-114.09  69.05\n";
	plotfile << "-113.89  70.12\n";
	plotfile << "-115.88  70.32\n";
	plotfile << "-116.10  71.32\n";
	plotfile << "-117.45  72.48\n";
	plotfile << "-113.53  72.44\n";
	plotfile << "-109.84  72.24\n";
	plotfile << "-106.62  71.71\n";
	plotfile << "-107.43  73.04 \n";
	plotfile << "\n";
	plotfile << "-120.96  74.29\n";
	plotfile << "-118.37  72.53\n";
	plotfile << "-123.06  71.18\n";
	plotfile << "-123.40  73.77\n";
	plotfile << "-120.93  74.27 \n";
	plotfile << "\n";
	plotfile << "-108.83  76.74\n";
	plotfile << "-106.25  75.54\n";
	plotfile << "-107.08  74.78\n";
	plotfile << "-112.99  74.16\n";
	plotfile << "-112.28  74.99\n";
	plotfile << "-116.04  75.33\n";
	plotfile << "-115.27  76.20\n";
	plotfile << "-110.95  75.56\n";
	plotfile << "-109.77  76.31\n";
	plotfile << "-108.82  76.70 \n";
	plotfile << "\n";
	plotfile << "-115.70  77.46\n";
	plotfile << "-118.10  76.30\n";
	plotfile << "-121.13  76.37\n";
	plotfile << "-116.04  77.28 \n";
	plotfile << "\n";
	plotfile << "-110.01  77.86\n";
	plotfile << "-112.36  77.68\n";
	plotfile << "-109.96  77.86 \n";
	plotfile << "\n";
	plotfile << "-109.60  78.48\n";
	plotfile << "-112.20  78.01\n";
	plotfile << "-109.60  78.48 \n";
	plotfile << "\n";
	plotfile << "-97.87  76.61\n";
	plotfile << "-99.21  75.31\n";
	plotfile << "-100.86  75.60\n";
	plotfile << "-99.40  76.26\n";
	plotfile << "-97.79  76.60 \n";
	plotfile << "\n";
	plotfile << "-94.72  75.53\n";
	plotfile << "-94.66  75.52 \n";
	plotfile << "\n";
	plotfile << "-104.10  79.01\n";
	plotfile << "-99.19  77.54\n";
	plotfile << "-103.22  78.08\n";
	plotfile << "-104.30  78.95 \n";
	plotfile << "\n";
	plotfile << "-93.74  77.52\n";
	plotfile << "-93.74  77.52 \n";
	plotfile << "\n";
	plotfile << "-96.88  78.50\n";
	plotfile << "-96.91  77.77\n";
	plotfile << "-96.94  78.48 \n";
	plotfile << "\n";
	plotfile << "-84.69  65.84\n";
	plotfile << "-81.58  63.87\n";
	plotfile << "-85.00  62.96\n";
	plotfile << "-84.63  65.71 \n";
	plotfile << "\n";
	plotfile << "-81.84  62.75\n";
	plotfile << "-82.01  62.63 \n";
	plotfile << "\n";
	plotfile << "-79.88  62.12\n";
	plotfile << "-79.88  62.12 \n";
	plotfile << "\n";
	plotfile << "-43.53  59.89\n";
	plotfile << "-45.29  60.67\n";
	plotfile << "-47.91  60.83\n";
	plotfile << "-49.90  62.41\n";
	plotfile << "-50.71  64.42\n";
	plotfile << "-51.39  64.94\n";
	plotfile << "-52.96  66.09\n";
	plotfile << "-53.62  67.19\n";
	plotfile << "-53.51  67.51\n";
	plotfile << "-51.84  68.65\n";
	plotfile << "-52.19  70.00\n";
	plotfile << "-51.85  71.03\n";
	plotfile << "-55.41  71.41\n";
	plotfile << "-54.63  72.97\n";
	plotfile << "-56.98  74.70\n";
	plotfile << "-61.95  76.09\n";
	plotfile << "-66.38  75.83\n";
	plotfile << "-71.13  77.00\n";
	plotfile << "-66.81  77.60\n";
	plotfile << "-70.78  77.78\n";
	plotfile << "-64.96  79.70\n";
	plotfile << "-63.38  81.16\n";
	plotfile << "-56.89  82.17\n";
	plotfile << "-48.18  82.15\n";
	plotfile << "-42.08  82.74\n";
	plotfile << "-38.02  83.54\n";
	plotfile << "-23.96  82.94\n";
	plotfile << "-25.97  81.97\n";
	plotfile << "-25.99  80.64\n";
	plotfile << "-13.57  80.97\n";
	plotfile << "-16.60  80.16\n";
	plotfile << "-19.82  78.82\n";
	plotfile << "-18.80  77.54\n";
	plotfile << "-21.98  76.46\n";
	plotfile << "-20.69  75.12\n";
	plotfile << "-21.78  74.40\n";
	plotfile << "-24.10  73.69\n";
	plotfile << "-26.54  73.08\n";
	plotfile << "-24.63  72.69\n";
	plotfile << "-21.84  71.69\n";
	plotfile << "-24.62  71.24\n";
	plotfile << "-27.16  70.89\n";
	plotfile << "-27.21  70.00\n";
	plotfile << "-24.10  69.35\n";
	plotfile << "-28.35  68.43\n";
	plotfile << "-32.48  68.56\n";
	plotfile << "-35.26  66.26\n";
	plotfile << "-37.90  65.90\n";
	plotfile << "-40.04  65.00\n";
	plotfile << "-40.49  64.04\n";
	plotfile << "-42.01  63.14\n";
	plotfile << "-42.88  61.15\n";
	plotfile << "-43.09  60.07\n";
	plotfile << "-43.56  59.90 \n";
	plotfile << "\n";
	plotfile << "-16.26  66.41\n";
	plotfile << "-15.32  64.29\n";
	plotfile << "-20.14  63.47\n";
	plotfile << "-21.76  64.21\n";
	plotfile << "-21.33  64.97\n";
	plotfile << "-23.04  65.62\n";
	plotfile << "-21.76  66.26\n";
	plotfile << "-18.77  66.12\n";
	plotfile << "-16.23  66.35 \n";
	plotfile << "\n";
	plotfile << "  0.56  51.47\n";
	plotfile << " -1.71  54.94\n";
	plotfile << " -3.41  57.52\n";
	plotfile << " -5.42  58.14\n";
	plotfile << " -5.77  55.59\n";
	plotfile << " -3.48  54.82\n";
	plotfile << " -4.68  52.88\n";
	plotfile << " -2.68  51.58\n";
	plotfile << " -3.80  50.08\n";
	plotfile << "  1.26  51.14\n";
	plotfile << "  0.65  51.41 \n";
	plotfile << "\n";
	plotfile << " -7.17  54.91\n";
	plotfile << " -9.97  53.47\n";
	plotfile << " -8.52  51.76\n";
	plotfile << " -5.69  54.79\n";
	plotfile << " -7.34  55.25 \n";
	plotfile << "\n";
	plotfile << " -1.33  60.66\n";
	plotfile << " -1.17  60.38 \n";
	plotfile << "\n";
	plotfile << " -6.18  58.44\n";
	plotfile << " -6.09  58.36 \n";
	plotfile << "\n";
	plotfile << " -6.47  57.58\n";
	plotfile << " -6.33  57.54 \n";
	plotfile << "\n";
	plotfile << " -7.30  57.54 \n";
	plotfile << "\n";
	plotfile << " -7.46  57.05 \n";
	plotfile << "\n";
	plotfile << " -6.54  56.94 \n";
	plotfile << "\n";
	plotfile << " -6.00  55.94 \n";
	plotfile << "\n";
	plotfile << " -5.09  55.55 \n";
	plotfile << "\n";
	plotfile << " -4.44  54.38\n";
	plotfile << " -4.30  54.19 \n";
	plotfile << "\n";
	plotfile << " -8.08  71.02\n";
	plotfile << " -8.21  70.86 \n";
	plotfile << "\n";
	plotfile << " 16.92  79.52\n";
	plotfile << " 22.26  78.46\n";
	plotfile << " 16.86  76.41\n";
	plotfile << " 16.00  77.39\n";
	plotfile << " 16.03  77.92\n";
	plotfile << " 16.81  79.50 \n";
	plotfile << "\n";
	plotfile << " 14.71  79.40\n";
	plotfile << " 16.05  79.12\n";
	plotfile << " 14.02  77.80\n";
	plotfile << " 13.56  78.46\n";
	plotfile << " 12.63  79.26\n";
	plotfile << " 14.68  79.40 \n";
	plotfile << "\n";
	plotfile << " 22.01  78.24\n";
	plotfile << " 21.86  78.23 \n";
	plotfile << "\n";
	plotfile << " 21.54  77.75\n";
	plotfile << " 23.88  77.26\n";
	plotfile << " 21.53  77.67\n";
	plotfile << " 22.79  77.79 \n";
	plotfile << "\n";
	plotfile << " 23.50  79.97\n";
	plotfile << " 28.24  79.54\n";
	plotfile << " 20.85  78.94\n";
	plotfile << " 19.00  79.34\n";
	plotfile << " 21.05  79.88\n";
	plotfile << " 23.41  79.96 \n";
	plotfile << "\n";
	plotfile << " 46.98  80.23\n";
	plotfile << " 43.13  79.97\n";
	plotfile << " 47.18  80.22\n";
	plotfile << "\n";
	plotfile << " 50.43  80.19\n";
	plotfile << " 50.55  79.88\n";
	plotfile << " 47.77  79.86\n";
	plotfile << " 50.45  80.14\n";
	plotfile << "\n";
	plotfile << " 61.79  80.18\n";
	plotfile << " 61.79  80.18\n";
	plotfile << "\n";
	plotfile << " 65.08  80.69\n";
	plotfile << " 64.27  80.59\n";
	plotfile << " 65.13  80.68\n";
	plotfile << "\n";
	plotfile << " -5.13  35.66\n";
	plotfile << "  4.06  36.63\n";
	plotfile << " 10.40  37.12\n";
	plotfile << " 11.36  33.61\n";
	plotfile << " 20.10  30.10\n";
	plotfile << " 23.49  32.17\n";
	plotfile << " 31.65  30.80\n";
	plotfile << " 35.76  23.74\n";
	plotfile << " 39.75  14.82\n";
	plotfile << " 42.93  11.34\n";
	plotfile << " 51.52  11.45\n";
	plotfile << " 49.82   6.99\n";
	plotfile << " 43.13  -0.62\n";
	plotfile << " 39.15  -7.58\n";
	plotfile << " 40.37 -13.20\n";
	plotfile << " 37.74 -18.17\n";
	plotfile << " 35.33 -22.71\n";
	plotfile << " 32.84 -28.15\n";
	plotfile << " 26.50 -34.39\n";
	plotfile << " 19.55 -35.51\n";
	plotfile << " 17.50 -30.88\n";
	plotfile << " 12.24 -18.75\n";
	plotfile << " 13.89 -12.81\n";
	plotfile << " 12.05  -5.55\n";
	plotfile << "  9.67   0.14\n";
	plotfile << "  7.19   3.79\n";
	plotfile << "  1.74   5.39\n";
	plotfile << " -4.77   4.59\n";
	plotfile << "-12.00   6.75\n";
	plotfile << "-15.54  10.98\n";
	plotfile << "-16.33  15.50\n";
	plotfile << "-16.10  22.29\n";
	plotfile << "-12.90  27.12\n";
	plotfile << " -9.52  31.09\n";
	plotfile << " -5.41  35.58\n";
	plotfile << "\n";
	plotfile << " 33.71   0.00\n";
	plotfile << " 33.48  -3.42\n";
	plotfile << " 33.34  -0.20\n";
	plotfile << " 33.71   0.00\n";
	plotfile << "\n";
	plotfile << " 49.30 -12.50\n";
	plotfile << " 49.28 -18.79\n";
	plotfile << " 43.95 -25.50\n";
	plotfile << " 44.37 -20.08\n";
	plotfile << " 46.34 -16.31\n";
	plotfile << " 47.91 -14.08\n";
	plotfile << " 49.30 -12.50\n";
	plotfile << "\n";
	plotfile << "178.88  69.10\n";
	plotfile << "181.20  68.42\n";
	plotfile << "183.52  67.78\n";
	plotfile << "188.87  66.38\n";
	plotfile << "186.54  64.74\n";
	plotfile << "182.87  65.63\n";
	plotfile << "180.13  65.14\n";
	plotfile << "179.48  64.88\n";
	plotfile << "178.20  64.29\n";
	plotfile << "177.46  62.62\n";
	plotfile << "170.42  60.17\n";
	plotfile << "164.48  59.89\n";
	plotfile << "162.92  57.34\n";
	plotfile << "161.82  54.88\n";
	plotfile << "156.42  51.09\n";
	plotfile << "156.40  57.76\n";
	plotfile << "163.79  61.73\n";
	plotfile << "159.90  60.73\n";
	plotfile << "156.81  61.68\n";
	plotfile << "153.83  59.10\n";
	plotfile << "148.57  59.46\n";
	plotfile << "140.77  58.39\n";
	plotfile << "137.10  54.07\n";
	plotfile << "140.72  52.43\n";
	plotfile << "138.77  47.30\n";
	plotfile << "129.92  42.04\n";
	plotfile << "128.33  38.46\n";
	plotfile << "126.15  35.18\n";
	plotfile << "125.12  39.08\n";
	plotfile << "121.62  40.15\n";
	plotfile << "117.58  38.21\n";
	plotfile << "121.77  36.90\n";
	plotfile << "120.73  32.65\n";
	plotfile << "121.28  30.25\n";
	plotfile << "118.83  24.93\n";
	plotfile << "112.69  21.81\n";
	plotfile << "108.53  21.73\n";
	plotfile << "107.55  16.34\n";
	plotfile << "107.32  10.45\n";
	plotfile << "104.39  10.37\n";
	plotfile << "100.01  13.52\n";
	plotfile << "100.26   8.30\n";
	plotfile << "103.22   1.56\n";
	plotfile << " 98.21   9.17\n";
	plotfile << " 97.66  15.36\n";
	plotfile << " 94.21  17.79\n";
	plotfile << " 90.05  21.74\n";
	plotfile << " 90.06  21.03\n";
	plotfile << " 82.06  15.95\n";
	plotfile << " 80.05  11.72\n";
	plotfile << " 76.41   8.60\n";
	plotfile << " 72.79  17.43\n";
	plotfile << " 72.02  20.00\n";
	plotfile << " 68.98  21.99\n";
	plotfile << " 64.62  24.41\n";
	plotfile << " 57.83  24.77\n";
	plotfile << " 53.11  26.20\n";
	plotfile << " 49.67  29.41\n";
	plotfile << " 50.96  25.15\n";
	plotfile << " 54.33  23.44\n";
	plotfile << " 59.03  22.57\n";
	plotfile << " 57.87  18.86\n";
	plotfile << " 52.95  15.74\n";
	plotfile << " 47.26  12.96\n";
	plotfile << " 42.75  14.68\n";
	plotfile << " 39.93  19.61\n";
	plotfile << " 36.92  25.78\n";
	plotfile << " 33.30  28.46\n";
	plotfile << " 32.60  30.63\n";
	plotfile << " 32.18  30.58\n";
	plotfile << " 36.08  35.03\n";
	plotfile << " 32.53  36.17\n";
	plotfile << " 27.77  36.94\n";
	plotfile << " 26.51  39.18\n";
	plotfile << " 31.54  40.82\n";
	plotfile << " 38.53  40.48\n";
	plotfile << " 40.35  43.17\n";
	plotfile << " 39.88  46.45\n";
	plotfile << " 35.18  44.99\n";
	plotfile << " 33.50  44.96\n";
	plotfile << " 30.24  45.14\n";
	plotfile << " 28.70  41.48\n";
	plotfile << " 26.55  39.84\n";
	plotfile << " 23.62  39.67\n";
	plotfile << " 23.80  37.34\n";
	plotfile << " 21.90  36.92\n";
	plotfile << " 18.79  42.02\n";
	plotfile << " 14.52  44.31\n";
	plotfile << " 14.58  42.25\n";
	plotfile << " 18.32  39.57\n";
	plotfile << " 16.05  39.35\n";
	plotfile << " 11.52  42.36\n";
	plotfile << "  6.87  43.08\n";
	plotfile << "  2.80  41.09\n";
	plotfile << " -1.11  37.14\n";
	plotfile << " -6.24  36.70\n";
	plotfile << " -8.67  39.57\n";
	plotfile << " -6.51  43.13\n";
	plotfile << " -0.84  45.55\n";
	plotfile << " -3.93  48.40\n";
	plotfile << "  0.48  49.09\n";
	plotfile << "  4.20  51.29\n";
	plotfile << "  6.44  52.92\n";
	plotfile << "  8.42  55.94\n";
	plotfile << " 11.72  55.49\n";
	plotfile << " 11.73  53.66\n";
	plotfile << " 16.78  54.14\n";
	plotfile << " 21.40  56.32\n";
	plotfile << " 24.67  57.20\n";
	plotfile << " 28.94  59.18\n";
	plotfile << " 24.16  59.52\n";
	plotfile << " 22.07  62.66\n";
	plotfile << " 23.76  65.35\n";
	plotfile << " 18.70  62.54\n";
	plotfile << " 19.11  59.67\n";
	plotfile << " 18.40  58.54\n";
	plotfile << " 15.34  55.73\n";
	plotfile << " 11.74  58.08\n";
	plotfile << "  8.37  57.68\n";
	plotfile << "  5.80  59.20\n";
	plotfile << "  7.38  60.86\n";
	plotfile << "  7.51  61.86\n";
	plotfile << "  9.62  62.99\n";
	plotfile << " 13.37  65.46\n";
	plotfile << " 15.46  67.12\n";
	plotfile << " 18.54  68.62\n";
	plotfile << " 22.32  69.64\n";
	plotfile << " 24.77  70.17\n";
	plotfile << " 25.93  69.79\n";
	plotfile << " 28.56  70.46\n";
	plotfile << " 29.75  69.76\n";
	plotfile << " 33.83  69.11\n";
	plotfile << " 41.90  66.85\n";
	plotfile << " 35.14  66.25\n";
	plotfile << " 33.30  66.07\n";
	plotfile << " 35.46  64.15\n";
	plotfile << " 37.68  64.03\n";
	plotfile << " 41.71  64.09\n";
	plotfile << " 44.80  65.58\n";
	plotfile << " 44.87  68.16\n";
	plotfile << " 45.92  66.83\n";
	plotfile << " 51.79  67.85\n";
	plotfile << " 53.70  67.89\n";
	plotfile << " 59.68  68.09\n";
	plotfile << " 65.07  69.08\n";
	plotfile << " 68.56  69.19\n";
	plotfile << " 68.38  70.97\n";
	plotfile << " 73.03  71.62\n";
	plotfile << " 73.80  68.29\n";
	plotfile << " 69.42  66.45\n";
	plotfile << " 73.43  66.36\n";
	plotfile << " 77.51  68.36\n";
	plotfile << " 80.74  66.74\n";
	plotfile << " 75.27  68.67\n";
	plotfile << " 75.11  71.80\n";
	plotfile << " 78.62  70.56\n";
	plotfile << " 78.43  71.90\n";
	plotfile << " 82.72  71.23\n";
	plotfile << " 84.25  70.03\n";
	plotfile << " 81.40  72.76\n";
	plotfile << " 86.50  74.01\n";
	plotfile << " 87.68  74.78\n";
	plotfile << " 90.25  75.23\n";
	plotfile << " 89.68  75.57\n";
	plotfile << " 95.12  75.95\n";
	plotfile << " 99.69  76.09\n";
	plotfile << "104.10  77.52\n";
	plotfile << "106.34  76.40\n";
	plotfile << "112.99  75.60\n";
	plotfile << "107.88  73.72\n";
	plotfile << "110.43  73.71\n";
	plotfile << "113.34  73.37\n";
	plotfile << "123.10  73.28\n";
	plotfile << "128.94  73.02\n";
	plotfile << "126.10  72.24\n";
	plotfile << "130.53  70.86\n";
	plotfile << "135.49  71.51\n";
	plotfile << "139.60  72.23\n";
	plotfile << "146.04  72.39\n";
	plotfile << "146.92  72.21\n";
	plotfile << "150.77  71.28\n";
	plotfile << "159.92  70.14\n";
	plotfile << "167.68  69.63\n";
	plotfile << "170.20  69.99\n";
	plotfile << "178.88  69.10\n";
	plotfile << "\n";
	plotfile << " 68.33  76.71\n";
	plotfile << " 66.03  75.62\n";
	plotfile << " 59.10  74.11\n";
	plotfile << " 54.92  73.03\n";
	plotfile << " 56.67  74.10\n";
	plotfile << " 58.56  75.09\n";
	plotfile << " 63.86  75.87\n";
	plotfile << " 68.19  76.70\n";
	plotfile << "\n";
	plotfile << " 53.04  72.57\n";
	plotfile << " 58.29  70.39\n";
	plotfile << " 55.03  70.78\n";
	plotfile << " 53.44  72.26\n";
	plotfile << " 53.63  72.61\n";
	plotfile << "\n";
	plotfile << " 52.22  46.50\n";
	plotfile << " 51.73  44.73\n";
	plotfile << " 52.56  41.80\n";
	plotfile << " 53.43  40.40\n";
	plotfile << " 54.22  37.86\n";
	plotfile << " 49.04  38.45\n";
	plotfile << " 48.17  42.76\n";
	plotfile << " 49.33  45.64\n";
	plotfile << " 52.22  46.50\n";
	plotfile << "\n";
	plotfile << " 62.32  46.32\n";
	plotfile << " 60.32  43.06\n";
	plotfile << " 59.57  45.58\n";
	plotfile << " 61.94  46.33\n";
	plotfile << "\n";
	plotfile << " 79.55  46.12\n";
	plotfile << " 74.30  44.44\n";
	plotfile << " 78.62  45.79\n";
	plotfile << " 79.66  46.07\n";
	plotfile << "\n";
	plotfile << " 76.81  41.96\n";
	plotfile << " 76.73  41.86\n";
	plotfile << "\n";
	plotfile << " 35.15  35.15\n";
	plotfile << " 34.61  34.84\n";
	plotfile << " 35.18  35.17\n";
	plotfile << "\n";
	plotfile << " 23.84  35.33\n";
	plotfile << " 24.30  34.91\n";
	plotfile << " 24.09  35.39\n";
	plotfile << "\n";
	plotfile << " 15.54  37.89\n";
	plotfile << " 13.47  37.89\n";
	plotfile << " 15.54  37.89\n";
	plotfile << "\n";
	plotfile << "  9.56  40.95\n";
	plotfile << "  8.46  39.99\n";
	plotfile << "  9.12  40.69\n";
	plotfile << "\n";
	plotfile << "  9.72  42.60\n";
	plotfile << "  9.54  42.35\n";
	plotfile << "\n";
	plotfile << " 80.60   8.95\n";
	plotfile << " 79.73   5.96\n";
	plotfile << " 80.10   8.30\n";
	plotfile << "\n";
	plotfile << " 11.04  57.44\n";
	plotfile << " 10.67  57.25\n";
	plotfile << "\n";
	plotfile << "-77.92  24.67\n";
	plotfile << "-77.98  24.22\n";
	plotfile << "\n";
	plotfile << "-77.61  23.62\n";
	plotfile << "-77.18  23.64\n";
	plotfile << "\n";
	plotfile << "-75.55  24.13\n";
	plotfile << "-75.41  24.31\n";
	plotfile << "\n";
	plotfile << "-91.40  -0.17\n";
	plotfile << "-91.52  -0.26\n";
	plotfile << "\n";
	plotfile << "-60.25  46.68\n";
	plotfile << "-60.71  46.33\n";
	plotfile << "\n";
	plotfile << "-63.89  49.47\n";
	plotfile << "-63.45  49.43\n";
	plotfile << "\n";
	plotfile << "142.53 -10.60\n";
	plotfile << "145.62 -16.34\n";
	plotfile << "149.79 -22.09\n";
	plotfile << "153.21 -26.82\n";
	plotfile << "150.52 -35.19\n";
	plotfile << "145.60 -38.53\n";
	plotfile << "140.13 -37.69\n";
	plotfile << "137.34 -34.77\n";
	plotfile << "135.76 -34.56\n";
	plotfile << "131.50 -31.34\n";
	plotfile << "121.72 -33.65\n";
	plotfile << "115.62 -33.25\n";
	plotfile << "114.09 -26.01\n";
	plotfile << "114.88 -21.27\n";
	plotfile << "122.34 -18.13\n";
	plotfile << "125.32 -14.53\n";
	plotfile << "128.39 -14.90\n";
	plotfile << "132.35 -11.42\n";
	plotfile << "136.16 -12.43\n";
	plotfile << "138.07 -16.45\n";
	plotfile << "142.25 -10.78\n";
	plotfile << "\n";
	plotfile << "144.72 -40.68\n";
	plotfile << "148.32 -42.14\n";
	plotfile << "145.57 -42.77\n";
	plotfile << "146.47 -41.19\n";
	plotfile << "\n";
	plotfile << "172.86 -34.23\n";
	plotfile << "176.10 -37.52\n";
	plotfile << "177.06 -39.49\n";
    plotfile << "174.77 -38.03\n";
	plotfile << "172.83 -34.27\n";
	plotfile << "\n";
	plotfile << "172.36 -40.53\n";
	plotfile << "172.92 -43.81\n";
	plotfile << "168.41 -46.13\n";
	plotfile << "170.26 -43.21\n";
	plotfile << "173.69 -40.94\n";
	plotfile << "\n";
	plotfile << "150.74 -10.18\n";
	plotfile << "143.04  -8.26\n";
	plotfile << "138.48  -6.97\n";
	plotfile << "131.95  -2.94\n";
	plotfile << "130.91  -1.35\n";
	plotfile << "134.38  -2.64\n";
	plotfile << "141.24  -2.62\n";
	plotfile << "148.19  -8.15\n";
	plotfile << "150.75 -10.27\n";
	plotfile << "\n";
	plotfile << "117.24   7.01\n";
	plotfile << "117.90   0.76\n";
	plotfile << "113.89  -3.50\n";
	plotfile << "109.44  -0.82\n";
	plotfile << "113.13   3.38\n";
	plotfile << "117.24   7.01\n";
	plotfile << "\n";
	plotfile << " 95.31   5.75\n";
	plotfile << "102.32   1.40\n";
	plotfile << "106.03  -2.98\n";
	plotfile << "101.46  -2.81\n";
	plotfile << " 95.20   5.73\n";
	plotfile << "\n";
	plotfile << "140.91  41.53\n";
	plotfile << "140.79  35.75\n";
	plotfile << "136.82  34.56\n";
	plotfile << "133.56  34.72\n";
	plotfile << "132.49  35.41\n";
	plotfile << "136.73  37.20\n";
	plotfile << "139.82  40.00\n";
	plotfile << "140.68  41.43\n";
	plotfile << "\n";
	plotfile << "133.71  34.30\n";
	plotfile << "131.41  31.58\n";
	plotfile << "129.38  33.10\n";
	plotfile << "133.90  34.37\n";
	plotfile << "\n";
	plotfile << "141.89  45.50\n";
	plotfile << "144.12  42.92\n";
	plotfile << "140.30  41.64\n";
	plotfile << "141.53  45.30\n";
	plotfile << "141.89  45.53\n";
	plotfile << "\n";
	plotfile << "142.57  54.36\n";
	plotfile << "143.64  49.19\n";
	plotfile << "141.99  45.88\n";
	plotfile << "141.92  50.85\n";
	plotfile << "142.60  54.34\n";
	plotfile << "\n";
	plotfile << "121.92  25.48\n";
	plotfile << "120.53  24.70\n";
	plotfile << "121.70  25.51\n";
	plotfile << "\n";
	plotfile << "110.81  20.07\n";
	plotfile << "109.20  19.66\n";
	plotfile << "110.81  20.07\n";
	plotfile << "\n";
	plotfile << "106.51  -6.16\n";
	plotfile << "114.15  -7.72\n";
	plotfile << "108.71  -7.89\n";
	plotfile << "106.51  -6.16\n";
	plotfile << "\n";
	plotfile << "164.27 -20.01\n";
	plotfile << "164.16 -20.27\n";
	plotfile << "\n";
	plotfile << "178.61 -17.04\n";
	plotfile << "178.61 -17.04\n";
	plotfile << "\n";
	plotfile << "179.45 -16.43\n";
	plotfile << "179.35 -16.43\n";
	plotfile << "\n";
	plotfile << "-172.55 -13.39\n";
	plotfile << "-172.61 -13.78\n";
	plotfile << "\n";
	plotfile << "122.26  18.67\n";
	plotfile << "123.05  13.86\n";
	plotfile << "120.73  13.80\n";
	plotfile << "120.43  16.43\n";
	plotfile << "121.72  18.40\n";
	plotfile << "\n";
	plotfile << "125.34   9.79\n";
	plotfile << "125.56   6.28\n";
	plotfile << "122.38   7.00\n";
	plotfile << "125.10   9.38\n";
	plotfile << "\n";
	plotfile << "119.64  11.35\n";
	plotfile << "118.81  10.16\n";
	plotfile << "119.59  10.86\n";
	plotfile << "119.64  11.35\n";
	plotfile << "\n";
	plotfile << "-179.87  65.14\n";
	plotfile << "-177.13  65.63\n";
	plotfile << "-173.46  64.74\n";
	plotfile << "-171.13  66.38\n";
	plotfile << "-176.48  67.78\n";
	plotfile << "-178.80  68.42\n";
	plotfile << "\n";
	plotfile << "101.96  79.08\n";
	plotfile << "101.31  77.86\n";
	plotfile << "101.22  79.04\n";
	plotfile << "\n";
	plotfile << " 94.29  79.29\n";
	plotfile << " 95.31  78.68\n";
	plotfile << "100.02  79.43\n";
	plotfile << " 97.26  79.62\n";
	plotfile << " 95.44  79.65\n";
	plotfile << "\n";
	plotfile << " 95.46  80.62\n";
	plotfile << " 92.39  79.66\n";
	plotfile << " 95.07  80.54\n";
	plotfile << "\n";
	plotfile << "138.54  76.05\n";
	plotfile << "144.93  75.45\n";
	plotfile << "140.30  74.99\n";
	plotfile << "137.27  75.44\n";
	plotfile << "138.29  75.98\n";
	plotfile << "\n";
	plotfile << "146.08  75.29\n";
	plotfile << "147.75  74.73\n";
	plotfile << "145.85  75.06\n";
	plotfile << "\n";
	plotfile << "141.44  73.88\n";
	plotfile << "141.48  73.84\n";
	plotfile << "\n";
	plotfile << "  0.01 -71.68\n";
	plotfile << "  6.57 -70.57\n";
	plotfile << " 15.04 -70.44\n";
	plotfile << " 25.10 -70.75\n";
	plotfile << " 33.37 -69.10\n";
	plotfile << " 38.46 -69.77\n";
	plotfile << " 42.85 -68.16\n";
	plotfile << " 46.59 -67.23\n";
	plotfile << " 49.35 -66.96\n";
	plotfile << " 52.90 -65.97\n";
	plotfile << " 58.46 -67.20\n";
	plotfile << " 63.60 -67.58\n";
	plotfile << " 70.63 -68.41\n";
	plotfile << " 69.24 -70.36\n";
	plotfile << " 76.20 -69.44\n";
	plotfile << " 88.08 -66.64\n";
	plotfile << " 94.98 -66.52\n";
	plotfile << "101.53 -66.09\n";
	plotfile << "111.31 -65.91\n";
	plotfile << "118.64 -66.87\n";
	plotfile << "126.24 -66.24\n";
	plotfile << "133.09 -66.18\n";
	plotfile << "139.85 -66.72\n";
	plotfile << "146.86 -67.96\n";
	plotfile << "153.65 -68.82\n";
	plotfile << "159.94 -69.57\n";
	plotfile << "164.10 -70.67\n";
	plotfile << "170.19 -71.94\n";
	plotfile << "165.68 -74.64\n";
	plotfile << "163.82 -77.60\n";
	plotfile << "162.10 -78.95\n";
	plotfile << "166.72 -82.84\n";
	plotfile << "175.58 -83.86\n";
	plotfile << "\n";
	plotfile << "-178.56 -84.37\n";
	plotfile << "-147.96 -85.40\n";
	plotfile << "-152.96 -81.12\n";
	plotfile << "-153.95 -79.50\n";
	plotfile << "-151.24 -77.48\n";
	plotfile << "-146.74 -76.44\n";
	plotfile << "-137.68 -75.16\n";
	plotfile << "-131.63 -74.63\n";
	plotfile << "-123.05 -74.41\n";
	plotfile << "-114.76 -73.97\n";
	plotfile << "-111.91 -75.41\n";
	plotfile << "-105.05 -74.77\n";
	plotfile << "-100.90 -74.21\n";
	plotfile << "-101.04 -73.18\n";
	plotfile << "-100.28 -73.06\n";
	plotfile << "-93.06 -73.33\n";
	plotfile << "-85.40 -73.18\n";
	plotfile << "-79.82 -73.04\n";
	plotfile << "-78.21 -72.52\n";
	plotfile << "-71.90 -73.41\n";
	plotfile << "-67.51 -71.10\n";
	plotfile << "-67.57 -68.92\n";
	plotfile << "-66.65 -66.83\n";
	plotfile << "-64.30 -65.28\n";
	plotfile << "-59.14 -63.74\n";
	plotfile << "-59.58 -64.37\n";
	plotfile << "-62.50 -65.94\n";
	plotfile << "-62.48 -66.66\n";
	plotfile << "-65.64 -68.02\n";
	plotfile << "-63.85 -69.07\n";
	plotfile << "-61.69 -70.87\n";
	plotfile << "-60.89 -72.71\n";
	plotfile << "-61.07 -74.30\n";
	plotfile << "-63.33 -75.88\n";
	plotfile << "-76.05 -77.06\n";
	plotfile << "-83.04 -77.12\n";
	plotfile << "-74.30 -80.83\n";
	plotfile << "-56.40 -82.14\n";
	plotfile << "-42.46 -81.65\n";
	plotfile << "-31.60 -80.17\n";
	plotfile << "-34.01 -79.20\n";
	plotfile << "-32.48 -77.28\n";
	plotfile << "-26.28 -76.18\n";
	plotfile << "-17.18 -73.45\n";
	plotfile << "-11.20 -72.01\n";
	plotfile << " -8.67 -71.98\n";
	plotfile << " -5.45 -71.45\n";
	plotfile << " -0.82 -71.74\n";
	plotfile << "  0.07 -71.68\n";
	plotfile << "\n";
	plotfile << "164.65 -77.89\n";
	plotfile << "170.95 -77.37\n";
	plotfile << "179.67 -78.25\n";
	plotfile << "\n";
	plotfile << "-178.74 -78.24\n";
	plotfile << "-165.76 -78.47\n";
	plotfile << "-158.42 -77.73\n";
	plotfile << "\n";
	plotfile << "-58.98 -64.63\n";
	plotfile << "-60.99 -68.62\n";
	plotfile << "-61.02 -71.70\n";
	plotfile << "\n";
	plotfile << "-62.01 -74.94\n";
	plotfile << "-52.00 -77.07\n";
	plotfile << "-42.23 -77.80\n";
	plotfile << "-36.22 -78.03\n";
	plotfile << "\n";
	plotfile << "-35.03 -77.81\n";
	plotfile << "-26.13 -75.54\n";
	plotfile << "-19.35 -73.04\n";
	plotfile << "-12.16 -71.86\n";
	plotfile << " -6.15 -70.65\n";
	plotfile << " -0.57 -69.14\n";
	plotfile << "  4.93 -70.25\n";
	plotfile << " 10.91 -69.99\n";
	plotfile << " 16.52 -69.87\n";
	plotfile << " 25.41 -70.22\n";
	plotfile << " 32.13 -69.29\n";
	plotfile << " 33.62 -69.58\n";
	plotfile << "\n";
	plotfile << " 70.56 -68.53\n";
	plotfile << " 73.91 -69.51\n";
	plotfile << "\n";
	plotfile << " 81.42 -67.87\n";
	plotfile << " 84.67 -66.41\n";
	plotfile << " 89.07 -66.73\n";
	plotfile << "\n";
	plotfile << "-135.79 -74.67\n";
	plotfile << "-124.34 -73.22\n";
	plotfile << "-116.65 -74.08\n";
	plotfile << "-109.93 -74.64\n";
	plotfile << "-105.36 -74.56\n";
	plotfile << "-105.83 -74.77\n";
	plotfile << "\n";
	plotfile << "-69.30 -70.06\n";
	plotfile << "-71.33 -72.68\n";
	plotfile << "-71.42 -71.85\n";
	plotfile << "-75.10 -71.46\n";
	plotfile << "-71.79 -70.55\n";
	plotfile << "-70.34 -69.26\n";
	plotfile << "-69.34 -70.13\n";
	plotfile << "\n";
	plotfile << "-49.20 -77.83\n";
	plotfile << "-44.59 -78.79\n";
	plotfile << "-44.14 -80.13\n";
	plotfile << "-59.04 -79.95\n";
	plotfile << "-49.28 -77.84\n";
	plotfile << "-48.24 -77.81\n";
	plotfile << "\n";
	plotfile << "-58.13 -80.12\n";
	plotfile << "-63.25 -80.20\n";
	plotfile << "-58.32 -80.12\n";
	plotfile << "\n";
	plotfile << "-163.64 -78.74\n";
	plotfile << "-161.20 -79.93\n";
	plotfile << "-163.62 -78.74\n";
	plotfile << "\n";
	plotfile << " 66.82  66.82\n";
	plotfile << " 66.82  66.82\n";
	plotfile.close();
}

int gnuplot::plot_world_map(
							string plot_script_filename,
							string plot_data_filename,
							string world_data_filename,
							string world_map_image_filename,
							string title,
							string subtitle, float indent, float vpos,
							vector<gridcell> &grid,
							int year,
							int reference_start_year,
							int reference_end_year,
							bool threed,
							float view_longitude, float view_latitude,
							bool show_cell_centres,
							int image_width,
							int image_height)
{
    ofstream plotfile;

	printf("Calculating reference temperatures...");

	for (int i = 0; i < (int)grid.size(); i++) {
		grid[i].update(reference_start_year,reference_end_year);
	}

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

	if (!ghcn::FileExists(world_data_filename)) {
		printf("done\nSaving world map data...");
		save_world_map(world_data_filename);
		printf("done");
	}

	printf("\nPlotting anomalies...");

    plotfile.open(plot_data_filename.c_str());
	for (int latitude=-89; latitude < 90; latitude++) {
		for (int longitude=-179; longitude < 180; longitude++) {

			int grid_cell_index =
				globalgrid::get_closest_grid_cell((float)-longitude, (float)latitude,grid);

			float anomaly =
				grid[grid_cell_index].get_anomaly(year,reference_start_year,reference_end_year);

			// 9999 indicates no data
			if (anomaly==-9999) anomaly=9999;
			
			if (show_cell_centres) {
				anomaly=0;
				for (int i = 0; i < (int)grid.size(); i++) {
					if (((int)(grid[i].latitude)==latitude) &&
						((int)(grid[i].longitude)==longitude)) {
						anomaly = 3;
						break;
					}
				}
			}

			plotfile << latitude << "  " << longitude << "      " << anomaly << "\n";
		}
		plotfile << "\n";
	}
	plotfile.close();


	string file_type = "png";
	if (world_map_image_filename != "") {
		if ((world_map_image_filename.substr((int)world_map_image_filename.size()-3,3) == "jpg") ||
			(world_map_image_filename.substr((int)world_map_image_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = world_map_image_filename.substr((int)world_map_image_filename.size()-3,3);
		}
	}

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";

	plotfile << "set output \"" << world_map_image_filename << "\"\n";
	plotfile << "set cbrange [-2:4]\n"; // Temperature range

	if (!threed) {
		// 2D map
		plotfile << "set border lw 1.5\n";
		plotfile << "set style line 1 lc rgb 'black' lt 1 lw 2\n";

		plotfile << "set rmargin screen 0.85\n";

		plotfile << "unset key\n";
		plotfile << "set tics scale 0.5\n";
		plotfile << "unset xtics\n";
		plotfile << "unset ytics\n";
		plotfile << "set xrange[-179:179]\n"; // longitude range
		plotfile << "set yrange[-89:89]\n"; // latitude range
		plotfile << "set format '%g'\n";
		plotfile << "set palette defined (0 \"blue\",17 \"#00ffff\",33 \"white\",50 \"yellow\",66 \"red\", 100 \"#990000\",101 \"grey\")\n";
		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";
		plotfile << "set title \"" << title << "\"\n";
		if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent - (0.34-0.17) << ", screen " << vpos << "\n";
		plotfile << "plot \"" << plot_data_filename << "\" u 2:1:3 w image, \"" << world_data_filename << "\" with lines linestyle 1\n";
	}
	else {
		// 3D map
		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_width << "\n";

		plotfile << "# color definitions\n";
		plotfile << "set border lw 1.5\n";
		plotfile << "set style line 1 lc rgb '#000000' lt 1 lw 2\n";
		plotfile << "set style line 2 lc rgb '#c0c0c0' lt 2 lw 1\n";
		plotfile << "unset key; unset border\n";
		plotfile << "set tics scale 0\n";
		plotfile << "set lmargin screen 0\n";
		plotfile << "set bmargin screen 0\n";
		plotfile << "set rmargin screen 1\n";
		plotfile << "set tmargin screen 1\n";
		plotfile << "set format ''\n";
		plotfile << "set mapping spherical\n";
		plotfile << "set angles degrees\n";
		plotfile << "set hidden3d front\n";
		plotfile << "# Set xy-plane to intersect z axis at -1 to avoid an offset between the lowest z\n";
		plotfile << "# value and the plane\n";
		plotfile << "set xyplane at -1\n";
		float lng = 90-view_longitude;
		if (lng<=0) lng += 360;
		if (lng>=360) lng -= 360;
		float lat = 90-view_latitude;
		if (lat<0) lat += 180;
		if (lat>=180) lat -= 180;
		plotfile << "set view " << lat << "," << lng << "\n";
		plotfile << "set parametric\n";
		plotfile << "set isosamples 25\n";
		plotfile << "set urange[0:360]\n";
		plotfile << "set vrange[-90:90]\n";
		plotfile << "set xrange[-1:1]\n";
		plotfile << "set yrange[-1:1]\n";
		plotfile << "set zrange[-1:1]\n";
		plotfile << "set palette defined (0 \"blue\",17 \"#00ffff\",33 \"white\",50 \"yellow\",66 \"red\",100 \"#990000\",101 \"grey\")\n";
		plotfile << "r = 0.99\n";
		plotfile << "set hidden3d front\n";
		plotfile << "splot \"" << plot_data_filename << "\" u 2:1:(1):3 w pm3d, r*cos(v)*cos(u),r*cos(v)*sin(u),r*sin(v) w l ls 2, \"" << world_data_filename << "\" u 1:2:(1) w l ls 1\n";
	}

	plotfile.close();

	printf("done\n");

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}

int gnuplot::plot_latitude_temperature(
									   string plot_script_filename,
									   string plot_data_filename,
									   string latitude_temperature_filename,
									   string title,
									   string subtitle, float indent, float vpos,
									   vector<tempdata> &data,
									   vector<stationdata> &stations,
									   int start_year,
									   int end_year,
									   string population_class,
									   bool show_minmax,
									   int image_width,
									   int image_height)
{
	int max_idx = (end_year - start_year+1);
	const int latitude_increment = 4;
	float latitude_temperature[180/latitude_increment];
	float latitude_temperature_min[180/latitude_increment];
	int latitude_temperature_min_hits[180/latitude_increment];
	float latitude_temperature_max[180/latitude_increment];
	int latitude_temperature_max_hits[180/latitude_increment];
	int latitude_temperature_hits[180/latitude_increment];
	float north=0;
	float west=0;
	float altitude=0;
	long long int prev_station=0;
	int idx=0;

	for (int i = 0; i < 180/latitude_increment; i++) {
		latitude_temperature[i]=0;
		latitude_temperature_hits[i]=0;

		latitude_temperature_min[i]=0;
		latitude_temperature_min_hits[i]=0;

		latitude_temperature_max[i]=0;
		latitude_temperature_max_hits[i]=0;
	}

	for (int i = 0; i < (int)data.size(); i++) {
		if ((data[i].year >= start_year) &&
			(data[i].year <= end_year)) {

			bool is_valid = true;

			if ((is_valid) && (population_class!="")) {
				is_valid = ghcn::station_has_property(
				    PROPERTY_POPULATION_CLASS,population_class,
					data[i].station, stations);
			}
			
			if ((is_valid) && (data[i].station!=prev_station)) {
				if (ghcn::get_station_location(
											   data[i].station,
											   stations,
											   north,west,altitude)) {
					idx =  ((int)(-north+0.5)+90)/latitude_increment;
					prev_station = data[i].station;
				}
				else {
					is_valid = false;
				}
			}

			if (is_valid) {
				latitude_temperature[idx] += data[i].temperature;
				latitude_temperature_hits[idx]++;
			}
		}
	}

	// find the minimum and maximum deviation
	if (show_minmax) {
		for (int i = 0; i < (int)data.size(); i++) {
			if ((data[i].year >= start_year) &&
				(data[i].year <= end_year)) {

				bool is_valid = true;

				if ((is_valid) && (population_class!="")) {
					is_valid = ghcn::station_has_property(
														  PROPERTY_POPULATION_CLASS,population_class,
														  data[i].station, stations);
				}

				if (data[i].station!=prev_station) {
					if (ghcn::get_station_location(
												   data[i].station,
												   stations,
												   north,west,altitude)) {
						idx =  ((int)(-north+0.5f)+90)/latitude_increment;
						prev_station = data[i].station;
					}
					else {
						is_valid = false;
					}
				}

				if (is_valid) {
					if (latitude_temperature_hits[idx]>0) {
						float av = latitude_temperature[idx]/latitude_temperature_hits[idx];
						if (data[i].temperature<av) {
							latitude_temperature_min[idx] += (av - data[i].temperature)*(av - data[i].temperature);
							latitude_temperature_min_hits[idx]++;
						}
						else {
							latitude_temperature_max[idx] += (data[i].temperature - av)*(data[i].temperature - av);
							latitude_temperature_max_hits[idx]++;
						}
					}
				}
			}
		}
	}

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    ofstream plotfile;
    plotfile.open(plot_data_filename.c_str());
	for (int i = 0; i < 180/latitude_increment; i++) {
		if (latitude_temperature_hits[i]>0) {
			float av = latitude_temperature[i]/latitude_temperature_hits[i];
			plotfile << (90-(i*latitude_increment)) << "    " << av;
			if (latitude_temperature_min_hits[i]>0) {
				plotfile << "    " << (av - sqrt((latitude_temperature_min[i]/latitude_temperature_min_hits[i])));
			}
			else {
				plotfile << "    0";
			}
			if (latitude_temperature_max_hits[i]>0) {
				plotfile << "    " << (av + sqrt((latitude_temperature_max[i]/latitude_temperature_max_hits[i])));
			}
			else {
				plotfile << "    0";
			}
			plotfile << "\n";
		}
    }
	plotfile.close();

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos << "\n";
	plotfile << "set yrange [-30:40]\n";
	plotfile << "set xrange [90:-90]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Latitude\"\n";
	plotfile << "set ylabel \"Average Temperature (Celcius)\"\n";
	plotfile << "set grid\n";
	plotfile << "set key right bottom\n";

	string file_type = "png";
	if (latitude_temperature_filename != "") {
		if ((latitude_temperature_filename.substr((int)latitude_temperature_filename.size()-3,3) == "jpg") ||
			(latitude_temperature_filename.substr((int)latitude_temperature_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = latitude_temperature_filename.substr((int)latitude_temperature_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";
		plotfile << "set output \"" << latitude_temperature_filename << "\"\n";
		if (!show_minmax) {
			plotfile << "plot \"" << plot_data_filename << "\" using 1:2 notitle with lines\n";
		}		
		else {
			plotfile << "plot \"" << plot_data_filename << "\" using 1:2 title \"Mean\" with lines";
			plotfile << ", \"" << plot_data_filename << "\" using 1:3 title \"Min\" with lines";
			plotfile << ", \"" << plot_data_filename << "\" using 1:4 title \"Max\" with lines\n";
		}
	}

	plotfile.close();

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}

int gnuplot::plot_altitude_temperature(
									   string plot_script_filename,
									   string plot_data_filename,
									   string altitude_temperature_filename,
									   string title,
									   string subtitle, float indent, float vpos,
									   vector<tempdata> &data,
									   vector<stationdata> &stations,
									   int start_year,
									   int end_year,
									   string population_class,
									   bool show_minmax,
									   int image_width,
									   int image_height)
{
	const int latitude_increment = 4;
	const int altitude_increment = 100;
	const int max_altitude = 4600;
	float latitude_temperature[180/latitude_increment];
	int latitude_temperature_hits[180/latitude_increment];
	float altitude_temperature[max_altitude/altitude_increment];
	int altitude_temperature_hits[max_altitude/altitude_increment];
	float north=0;
	float west=0;
	float altitude=0;
	long long int prev_station=0;
	int idx=0;
	float min_temperature = -1;
	float max_temperature = 1;

	for (int i = 0; i < 180/latitude_increment; i++) {
		latitude_temperature[i]=0;
		latitude_temperature_hits[i]=0;
	}
	for (int i = 0; i < max_altitude/altitude_increment; i++) {
		altitude_temperature[i]=0;
		altitude_temperature_hits[i]=0;
	}

	for (int i = 0; i < (int)data.size(); i++) {
		if ((data[i].year >= start_year) &&
			(data[i].year <= end_year)) {

			bool is_valid = true;

			if ((is_valid) && (population_class!="")) {
				is_valid = ghcn::station_has_property(
				    PROPERTY_POPULATION_CLASS,population_class,
					data[i].station, stations);
			}
			
			if ((is_valid) && (data[i].station!=prev_station)) {
				if (ghcn::get_station_location(
											   data[i].station,
											   stations,
											   north,west,altitude)) {
					idx =  ((int)(-north+0.5)+90)/latitude_increment;
					prev_station = data[i].station;
				}
				else {
					is_valid = false;
				}
			}

			if (is_valid) {
				latitude_temperature[idx] += data[i].temperature;
				latitude_temperature_hits[idx]++;
			}
		}
	}

	// latitude adjusted temperatures
	for (int i = 0; i < (int)data.size(); i++) {
		if ((data[i].year >= start_year) &&
			(data[i].year <= end_year)) {

			bool is_valid = true;

			if ((is_valid) && (population_class!="")) {
				is_valid = ghcn::station_has_property(
													  PROPERTY_POPULATION_CLASS,population_class,
													  data[i].station, stations);
			}

			if (data[i].station!=prev_station) {
				if (ghcn::get_station_location(
											   data[i].station,
											   stations,
											   north,west,altitude)) {
					idx =  ((int)(-north+0.5f)+90)/latitude_increment;
					prev_station = data[i].station;
				}
				else {
					is_valid = false;
				}
			}

			if (is_valid) {
				if ((latitude_temperature_hits[idx]>0) && (altitude<max_altitude) && (altitude>0)) {
					float alt_temp = data[i].temperature - (latitude_temperature[idx]/latitude_temperature_hits[idx]);
					int alt_idx = (int)(altitude/altitude_increment);
					altitude_temperature[alt_idx] += data[i].temperature - (latitude_temperature[idx]/latitude_temperature_hits[idx]);
					altitude_temperature_hits[alt_idx]++;
				}
			}
		}
	}

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    ofstream plotfile;
    plotfile.open(plot_data_filename.c_str());
	for (int i = 0; i < max_altitude/altitude_increment; i++) {
		if (altitude_temperature_hits[i]>0) {
			float av = altitude_temperature[i]/altitude_temperature_hits[i];
			if (av > max_temperature) max_temperature = av;
			if (av < min_temperature) min_temperature = av;
			plotfile << (i*altitude_increment) << "    " << av << "\n";
		}
    }
	plotfile.close();

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos << "\n";
	plotfile << "set yrange [" << min_temperature << ":" << max_temperature << "]\n";
	plotfile << "set xrange [0:" << max_altitude << "]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Altitude (Metres above Sea Level)\"\n";
	plotfile << "set ylabel \"Average Latitide Adjusted Temperature (Celcius)\"\n";
	plotfile << "set grid\n";
	plotfile << "set key right bottom\n";

	string file_type = "png";
	if (altitude_temperature_filename != "") {
		if ((altitude_temperature_filename.substr((int)altitude_temperature_filename.size()-3,3) == "jpg") ||
			(altitude_temperature_filename.substr((int)altitude_temperature_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = altitude_temperature_filename.substr((int)altitude_temperature_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";
		plotfile << "set output \"" << altitude_temperature_filename << "\"\n";
		if (!show_minmax) {
			plotfile << "plot \"" << plot_data_filename << "\" using 1:2 notitle with lines\n";
		}		
		else {
			plotfile << "plot \"" << plot_data_filename << "\" using 1:2 title \"Mean\" with lines";
			plotfile << ", \"" << plot_data_filename << "\" using 1:3 title \"Min\" with lines";
			plotfile << ", \"" << plot_data_filename << "\" using 1:4 title \"Max\" with lines\n";
		}
	}

	plotfile.close();

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}

int gnuplot::plot_station_altitudes(
								  string plot_script_filename,
								  string plot_data_filename,
								  string station_altitudes_filename,
								  string title,
								  string subtitle, float indent, float vpos,
								  vector<tempdata> &data,
								  vector<long long int> &country_codes,
								  vector<stationdata> &stations,
								  int start_year,
								  int end_year,
								  float min_temp, float max_temp,
								  string population_class,
								  vector<float> &area,
								  int image_width,
								  int image_height)
{
	const int max_altitude = 2500;
	const int altitude_increment = 50;
	int altitude_histogram[max_altitude/altitude_increment];
	long long int prev_station=0;
	float north=0,west=0,altitude=0;
	int idx=-1;

	for (int i = 0; i < max_altitude/altitude_increment; i++) {
		altitude_histogram[i]=0;
	}

	for (int i = 0; i < (int)data.size(); i++) {
		if ((data[i].year >= start_year) &&
			(data[i].year <= end_year) &&
			(data[i].temperature>=min_temp) &&
			(data[i].temperature<=max_temp)) {

			bool is_valid = true;
			if ((int)country_codes.size() > 0) {
                is_valid = ghcn::vector_contains(country_codes, (long long int)data[i].country_code);
			}

			if (is_valid) {
				if (data[i].station!=prev_station) {
					if (ghcn::get_station_location(
												   data[i].station,
												   stations,
												   north,west,altitude)) {
						if ((altitude>=0) && (altitude<max_altitude)) {
							idx = (int)(altitude/altitude_increment);

							if (area.size()==4) {
								if (!ghcn::station_within_area(data[i].station,stations,
															   area[0],area[1],area[2],area[3])) {
									idx=-1;
								}
							}

							if (population_class!="") {
								if (!ghcn::station_has_property(PROPERTY_POPULATION_CLASS,population_class,
																data[i].station, stations)) {
									idx=-1;
								}
							}

						}
						else {
							idx=-1;
						}
						prev_station = data[i].station;
					}
				}

				if (idx>-1) {
					altitude_histogram[idx]++;
				}

			}
		}
	}


	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    ofstream plotfile;
	int max=1;
    plotfile.open(plot_data_filename.c_str());
	for (int i = 0; i < max_altitude/altitude_increment; i++) {
		if (altitude_histogram[i]>max) max = altitude_histogram[i];
    	plotfile << (i*altitude_increment) << "    " << altitude_histogram[i]/1000.0f << "\n";
    }
	plotfile.close();

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos << "\n";
	plotfile << "set yrange [" << 0 << ":" << (max*105/100000) << "]\n";
	plotfile << "set xrange [0:" << max_altitude << "]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Altitude (Metres above Sea Level)\"\n";
	plotfile << "set ylabel \"Number of samples x 1000\"\n";
	plotfile << "set grid\n";
	plotfile << "set key right bottom\n";

	string file_type = "png";
	if (station_altitudes_filename != "") {
		if ((station_altitudes_filename.substr((int)station_altitudes_filename.size()-3,3) == "jpg") ||
			(station_altitudes_filename.substr((int)station_altitudes_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = station_altitudes_filename.substr((int)station_altitudes_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";
		plotfile << "set output \"" << station_altitudes_filename << "\"\n";
		plotfile << "plot " << "\"" << plot_data_filename << "\" using 1:2 notitle with lines\n";
	}

	plotfile.close();

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}

int gnuplot::plot_average_altitudes(
							string plot_script_filename,
							string plot_data_filename,
							string altitudes_filename,
							string title,
							string subtitle, float indent, float vpos,
							vector<tempdata> &data,
							vector<long long int> &country_codes,
							vector<stationdata> &stations,
							int start_year,
							int end_year,
							float min_temp, float max_temp,
							string population_class,
							vector<float> &area,
							bool show_running_average,
							int image_width,
							int image_height)
{
	const int max_altitude = 4600;
	int max_idx = (end_year - start_year+1);
	float altitude[max_idx];
	int altitude_hits[max_idx];
	long long int prev_station=0;
	float north=0,west=0,alt=0;
	int idx=-1;
	bool found = false;

	for (int i = 0; i < max_idx; i++) {
		altitude[i]=0;
		altitude_hits[i]=0;
	}

	for (int i = 0; i < (int)data.size(); i++) {
		if ((data[i].year >= start_year) &&
			(data[i].year <= end_year) &&
			(data[i].temperature>=min_temp) &&
			(data[i].temperature<=max_temp)) {

			bool is_valid = true;
			if ((int)country_codes.size() > 0) {
                is_valid = ghcn::vector_contains(country_codes, (long long int)data[i].country_code);
			}

			if (is_valid) {
				if (data[i].station!=prev_station) {
					if (ghcn::get_station_location(
												   data[i].station,
												   stations,
												   north,west,alt)) {
						if ((alt>=0) && (alt<max_altitude)) {
							found = true;

							if (area.size()==4) {
								if (!ghcn::station_within_area(data[i].station,stations,
															   area[0],area[1],area[2],area[3])) {
									found = false;
								}
							}

							if (population_class!="") {
								if (!ghcn::station_has_property(PROPERTY_POPULATION_CLASS,population_class,
																data[i].station, stations)) {
									found = false;
								}
							}

						}
						else {
							found = false;
						}
						prev_station = data[i].station;
					}
				}

				if (found) {
					idx = data[i].year - start_year;
					altitude[idx]+=alt;
					altitude_hits[idx]++;
				}

			}
		}
	}

	std::remove(plot_script_filename.c_str());
	std::remove(plot_data_filename.c_str());

    ofstream plotfile;
	int max=1;
	float running_average=-1;
    plotfile.open(plot_data_filename.c_str());
	for (int i = 0; i < max_idx; i++) {
		if (altitude_hits[i]>0) {
			altitude[i] /= altitude_hits[i];

			if (running_average<0) {
				running_average = altitude[i];
			}
			else {
				running_average = (running_average*0.8f) + (altitude[i]*0.2f);
			}

			if (show_running_average) {
				if (running_average > max) max = running_average;
				plotfile << (start_year + i) << "    " << running_average << "\n";
			}
			else {
				if (altitude[i] > max) max = altitude[i];
				plotfile << (start_year + i) << "    " << altitude[i] << "\n";
			}
		}
    }
	plotfile.close();

    plotfile.open(plot_script_filename.c_str());
	plotfile << "reset\n";
	plotfile << "set title \"" << title << "\"\n";
	if (subtitle != "") plotfile << "set label \"" << subtitle << "\" at screen " << indent << ", screen " << vpos << "\n";
	plotfile << "set yrange [0:" << (max*105/100) << "]\n";
	plotfile << "set xrange [" << start_year << ":" << end_year << "]\n";
	plotfile << "set lmargin 9\n";
	plotfile << "set rmargin 2\n";
	plotfile << "set xlabel \"Year\"\n";
	plotfile << "set ylabel \"Average Altitude (Metres above Sea Level)\"\n";
	plotfile << "set grid\n";
	plotfile << "set key right bottom\n";

	string file_type = "png";
	if (altitudes_filename != "") {
		if ((altitudes_filename.substr((int)altitudes_filename.size()-3,3) == "jpg") ||
			(altitudes_filename.substr((int)altitudes_filename.size()-4,4) == "jpeg")) {
			file_type = "jpeg";
		}
		else {
			file_type = altitudes_filename.substr((int)altitudes_filename.size()-3,3);
		}

		plotfile << "set terminal " << file_type << " size " << image_width << "," << image_height << "\n";
		plotfile << "set output \"" << altitudes_filename << "\"\n";
		plotfile << "plot " << "\"" << plot_data_filename << "\" using 1:2 notitle with lines\n";
	}

	plotfile.close();

	string command_str = "gnuplot " + plot_script_filename;
	return system(command_str.c_str());
}
